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This  research  has  dealt  with  optimal  control  and  adaptive 
parameter  identification  and  control  of  distributed  systems,  with 
primary  application  to  large  flexible  space  structures.  The 
research  has  developed  both  mathematical  theory  and  numerical 
methods  for  design  of  control  laws  and  compensators  for  complex 
distributed  systems  like  large  space  platforms  and  antennas. 
Approximation  theory  that  serves  as  the  basis  for  computer  aided 
design  of  control  systems  involving  partial  and  functional 
differential  equations  has  been  an  especially  improtant  part  of 
this  research  project.  Also,  variable-order  adaptive  parameter 
estimators  and  adaptive  controllers  have  been  developed  for 
distributed  systems,  and  the  parameter  estimators  (lattice 
filters)  have  been  applied  to  both  numerical  simulations  and 
experimental  structures. 


In  optimal  control  of  distributed  systems,  we  have  studied 
control  problems  for  infinite  dimensional  systems  with  quadratic 
performance  criteria.  These  problems  require  the  solution  of 
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infinite  dimensional  Riccati  operator  equations,  which  must  be 
solved  numerically  by  large-order  approximations.  Approximation 
theory  and  numerical  results  for  continuous-time  optimal  control 
are  given  in  [1,2,16]. 

We  have  collaborated  with  Professor  I.  Gary  Rosen  to  develop 
approximation  theory  and  numerical  methods  for  discrete-time 
infinite  dimensional  optimal  regulator  problems  and  to 
demonstrate  the  practical  implications  of  the  theory  and  the 
usefulness  of  the  numerical  methods  on  several  examples, 
involving  flexible  structures,  diffusion  processes  and  time-delay 
systems  [3,4,5,6,15].  In  these  numerical  examples,  we  have 
discovered  important  relationships  between  the  length  of  the 
sampling  interval  (the  time  step)  and  the  numerical  condition  of 
the  approxomating  finite-dimensional  discrete-time  control 
problems . 


In  connection  with  optimal  control  of  flexible  space 
structures,  we  have  investigated  integrated  control-structure 
design,  by  which  we  mean  simultaneous  optimal  design  of  structure 
and  controller.  We  have  combined  structural  weight  and 
robustness  with  respect  to  modelling  errors  in  an  overall 
objective  functional  for  integrated  control -structure 
optimization.  This  results  in  a  nonlinear  programming  problem 


whose  solution  we  have  demonstrated  in  several  numerical  examples 


[7,8,16,17]. 
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Another  area  of  primary  interest  in  this  project  has  been 
adaptive  parameter  identification  and  adaptive  control  of 
distributed  systems.  We  have  studied  the  application  of  modern 
least-square  identification  schemes  to  infinite  dimensional 
systems,  and  our  most  important  results  concern  lattice  filters 
for  estimating  coefficients  in  infinite  dimensional 
autoregressive-moving  average  (ARMA)  models  [9,10,11,12,15]. 

Because  lattice  filters  are  recursive  in  order  as  well  as  in 
time,  they  permit  the  order  of  the  model  to  vary  adaptively  as 
either  more  structural  modes  are  excited  or  as  higher  frequency 
modes  are  damped  out  while  lower  modes  remain  excited.  We  have 
developed  a  vector-channel  lattice  filter  [10,15],  which  is  an 
extension  of  previous  lattice  filters  and  which  is  important  for 
applications  with  several  measurements  from  the  same  structure. 
With  experimental  data  from  a  large  and  very  flexible  grid  at  the 
NASA  Langley  Research  Center,  we  have  demonstrated  the 
effectiveness  of  lattice  filters  for  adaptive  identification  and 
prediction  of  flexible  structures  [11]. 

We  have  begun  developing  approximation  theory  that  will 
predict  the  performance  of  lattice  filters  with  structures  having 
various  types  of  damping  [12].  We  expect  such  theory  to  be 
important  in  both  adaptive  identification  and  adaptive  control  of 
large  space  structures. 
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Also,  we  have  incorporated  our  adaptive  parameter  estimators 
into  variable-order  adaptive  control  algorithms  for  flexible 
space  structrues.  Our  main  application  so  far  has  been  to 
robotic  manipulators  with  flexible  links  [13,14,18,19].  We  plan 
to  develop  similar  control  algorithms  for  large  flexible 
platforms  and  antennas  with  structure  similar  to  the  flexible 
grid  on  which  [11]  is  based. 
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ABSTRACT 


We  consider  the  approximation  of  optimal  discrete-time  linear  quadratic  Gaussian  (LQG) 
compensators  for  distributed  parameter  control  systems  with  boundary  input  and  unbounded 
measurement.  Our  approach  applies  to  problems  that  can  be  formulated  in  a  state  space  on  which 
both  the  discrete-time  input  and  output  operators  are  continuous.  Approximating  compensators  are 
obtained  via  application  of  the  LQG  theory  and  associated  approximation  results  for  infinite 
dimensional  discrete-time  control  systems  with  bounded  input  and  output.  Numerical  results  for 
spline  and  modal  approximation  schemes  used  to  compute  optimal  compensators  for  a  one 
dimensional  heat  equation  with  either  Neumann  or  Dirichlet  boundary  control  and  pointwise 
measurement  of  temperature  are  presented  and  discussed. 


In  this  paper  we  develop  an  approximation  theory  for  the  computation  of  optimal  discrete-time 
linear  quadratic  Gaussian  (LQG)  compensators  (combined  feedback  control  law  and  state  estimator) 
for  distributed  parameter  systems  with  boundary  input  or  control  and  unbounded  output  or 
measurement.  In  a  continuous  time  setting,  boundary  input  typically  results  in  an  unbounded  input 
operator.  That  is,  the  system's  input  operator  maps  the  control  input  into  a  space  larger  than  the 
state  space  in  which  the  open-loop  system  is  usually  formulated.  In  the  discrete-time  case,  on  the 
other  hand,  for  a  wide  class  of  distributed  systems,  the  resulting  input  operator  is  bounded  on  the 
usual  underlying  state  space.  An  unbounded  output,  or  measurement,  operator  has  domain  smaller 
than  the  usual  open-loop  state  space. 

For  continuous  time  systems,  Pritchard  and  Salamon  (1987)  have  established  an  abstract 
semigroup  theoretic  framework  for  treating  the  linear  quadratic  regulator  problem  (control  only)  for 
infinite  dimensional  systems  with  unbounded  input  and  output  operators.  Their  approach  is  based 
upon  a  weak  or  distributional  formulation  of  the  Riccati  equations  which  characterize  the  optimal 
feedback  control  laws  in  an  appropriate  dual  space  .  Curtain  (1984)  provides  a  procedure  for  the 
design  of  finite  dimensional  compensators  for  parabolic  systems  with  unbounded  control  and 
observation.  In  (Curtain  and  Salamon,  1986)  a  finite  dimensional  compensator  design  procedure 
for  a  wider  class  of  infinite  dimensional  systems  with  unbounded  input  (but  bounded  output) 
including  hereditary  systems  with  control  delays  and  partial  differential  systems  with  boundary 
control  is  developed.  Lasiecka  and  Triggiani  have  looked  at  linear  regulator  problems  for  parabolic 
(1983a,  1987a)  and  hyperbolic  (1983b,  1986)  systems  with  boundary  control  and  obtained, 
among  other  things,  global  and  local  regularity  results  for  the  optimal  controls  and  state 
trajectories.  In  (Lasiecka  and  Triggiani,  1987b)  Galerkin  approximations  and  an  associated 
convergence  theory  for  closed-loop  solutions  to  regulator  problems  for  parabolic  systems  with 
Dirichlet  boundary  input  are  studied.  A  more  complete  survey  of  the  boundary  control  literature 
including  references  to  some  of  the  poineering  work  in  this  area  can  be  found  in  (Pritchard  and 
Salamon,  1987). 


In  our  treatment  here,  we  consider  the  discrete-time  problem  (i.e.  piecewise  constant  input  and 
sampled  output).  Our  interest  in  the  discrete-time  or  digital  formulation  is  motivated  by  1)  the  fact 
that  it  represents  a  more  accurate  or  realistic  description  of  how  the  linear-quadratic  theory  for 
distributed  systems  would  actually  be  applied  in  practice,  and  by  2)  how  the  boundedness  of  the 
discrete-time  input  operator  in  the  usual  underlying  state  space  facilitates  the  development  of  an 
approximation  theory  which  can  simultaneously  handle  both  unbounded  input  and  unbounded 
output.  Our  approach  is  based  upon  an  application  of  the  theory  we  developed  earlier  in  (Gibson 
and  Rosen,  1985  and  1986)  for  the  approximation  of  optimal  discrete-time  LQG  compensators  for 
infinite  dimensional  systems  with  bounded  input  and  output.  Our  results  are  applicable  to 
boundary  control  systems  in  which  a  restriction  of  the  state  transition  operator  and  the  discrete-time 
input  operator  are  bounded  on  a  space  on  which  the  output  operator  is  bounded  as  well.  To 
illustrate  our  approach,  in  this  paper  we  describe  in  detail  the  application  of  our  theory  to  a  one 
dimensional  heat  equation  with  either  Neumann  or  Dirichlet  boundary  control  and  pointwise 
measurement  of  temperature.  Elsewhere  (see  Gibson  and  Rosen,  1987)  we  have  applied  our 
results  to  develop  approximation  schemes  for  the  computation  of  optimal  LQG  compensators  for 
flexible  structures  (i.e.  Euler-Bemoulli  beams)  with  shear  force  input  at  the  boundary  and  a 
pointwise  measurement  of  strain. 

An  outline  of  the  remainder  of  the  paper  is  as  follows.  In  Section  2  we  describe  an  abstract 
framework  for  the  study  of  boundary  control  systems  and  their  discrete-time  formulation.  In 
Section  3  we  review  the  LQG  theory  for  infinite  dimensional  discrete-time  systems  and  associated 
abstract  approximation  results.  In  the  fourth  section,  we  discuss  spline  and  modal  subspace  based 
approximation  schemes  for  the  heat  equation  example.  Section  5  contains  some  concluding 
remarks. 


We  employ  a  semigroup  theoretic  formulation  that  has  been  used  previously  for  a  class  of 
abstract  boundary  control  systems.  See,  for  example,  (Curtain  and  Salamon,  1986).  Let  W,V  and 
H  be  Hilbert  spaces  with  W  and  V  densely  and  continuously  embedded  in  H.  We 
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consider  boundary  control  systems  of  the  form 


(2.1) 

✓—s 

II 

> 

t>0 

(2.2) 

w(0)  =  w0 

(2.3) 

Tw(t)  =  v(t), 

t>0 

(2.4) 

y(t)  =  Cw(t), 

t>0 

where  A  e  £(W,H),  the  boundary  input  operator  T  is  an  element  in  £(W,Rm)  and  the  output 
operator  C  is  an  element  in  £(V,RP).  Note  that  the  operator  A  need  not  be  the  Laplacian.  Our 
choice  of  A  to  denote  a  general,  most  often  differential,  operator  satisfying  the  conditions  set  forth 
below  is  consistent  with  the  notation  used  in  earlier  treatments  of  boundary  control  systems 
elsewhere  in  the  literature. 

We  assume  that  1)  T  is  surjective  and  its  null  space,  7l(T)  =  {9  e  W:  Hp  =  0},  is  dense  in 
H,  2)  the  operator  0,  defined  to  be  the  restriction  of  the  operator  A  to  7UD.  is  a  closed  operator 
on  H  and  has  non-empty  resolvent  set  and  3)  for  each  T  >  0,  all  wQ  e  W,  and  v  e  C1  (0,T;  Rm) 
with  rw0  =  v(0),  there  exists  a  unique  w  e  C([0,T];  W)  D  C'dO.TJ;  H)  which  depends 
continuously  on  wQ  and  v  and  which  satisfies  (2.1)-  (2.3)  for  each  t  e  [0,T].  It  then  follows  (see 
Hille  and  Phillips,  1957)  that  the  operator  Cf :  Dom  (Cl)  c  H  ->  H  given  by  Cup  =  Atp  for  cp  e 
Dom(C0  =  TUOis  the  infinitesimal  generator  of  a  £0  semigroup,  {ST(t):  t>  0),  of  bounded 
linear  operators  on  H. 

Define  the  space  Z  as  the  dual  of  Dom(Cf  )  where  the  norm  on  Dom  (Ct  )  is  taken  to  be  the 
graph  Hilbert  space  norm  associated  with  the  operator  Ct  .  Then  H  is  densely  and  continuously 
embedded  in  Z  and  { <T(t) :  t  >  0}  can  be  uniquely  extended  to  a  G0  semigroup  of  bounded  linear 
operators  on  Z.  Its  generator  is  the  extension  of  the  operator  Q  to  an  the  operator  C(  in  C(H,Z) 
given  by  (  Cf  <p)(\p)  =  <(p,  Cf*y>H  for  tp  e  H  and  \pe  Dom(Cf*). 

Since  T  was  assumed  to  be  a  surjection,  it  has  a  right  inverse.  Let  P :  Rm  — »  W  be  any  right 
inverse  of  T .  Since  Dom  (  P)  =  Rm,  we  have  P  e  £  (Rm,  W).  For  v  e  Rm,  we  define 
fa  6  Z(Rm,  Z)  by  fav  =  (A  -  Cf)  P  v.  If  P  and  P  are  two  distinct  right  inverses  of  T  then 
fa,(  P  -  P )  c  71(0-  Since  Ct  coincides  with  A  on  71(f) ,  it  follows  that  the  operator  fa  is 
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well  defined.  It  can  be  shown  (see  Curtain  and  Salamon,  1986)  that  for  each  w0e  H  and  v  E 
1^(0, T;  Rm)  there  exists  a  unique  w  e  C([0,T];  H)  fl  H1  (0,T;  Z)  which  depends  continuously  on 
w0  and  v  and  which  satisfies 

w(t)  =  Ciw(t)  +  Bv(t),  t  >  0 
w(0)  =  w0 

in  Z.  The  function  w  is  given  by 

t 

(2.5)  w(t)  =tT(t)w0  +  J  $r(t -  s)  Bv(s)ds,  t>0 

o 

and  is  referred  to  as  a  weak  solution  to  the  boundary  control  system  (2. 1)  -  (2.3). 

The  discrete-time  formulation  of  (2.1)  -  (2.4)  is  found  by  considering  piecewise-constant 
controls  of  the  form 

(2.6)  v(t)  =  uk,  t  e  [kx,  (k  +  l)x),  k  =  0,1,2,... 

where  x  denotes  the  length  of  the  sampling  interval.  Let  wk  =  w(kx) ,  k=  0,1,2,...  where  w(  ) 
is  the  unique  weak  solution  to  (2.1)  -  (2.3)  given  by  (2.5)  corresponding  to  w0  £  H  and  input  v 
given  by  (2.6).  (We  note  that  with  piecewise  constant  input  of  the  form  (2.6),  the  solution  w  is  in 
fact  a  strong  solution  on  each  subinterval  [kx,  (k+l)x].)  For  each  k  =  0,1,2,...  define 
zk  £  C([kx,  (k  +  l)xj;  H)  by  z^t)  =  w(t)  -  r^,  1  e  lkx>  <k  +  W  171611 

zk(t)  =  w(t)  =  Clw(t)  +  Buk 

=  ci  zk(t)  +  (Ci  +  Bnr^X 
=  Clz^t)  +  A  r*uk ,  t  £  (kx,  (k  +  l)x], 
z^k  x)  =  wk  -  ruk  . 


Therefore 


Wk+1  =  z^Cfk+Ox)  +  Puk 

=  <r (x)(wk  -  r-uk)  +  o' (s)  a  r+^ds  +  rfuk 

=  cr (x)wk  +  a  -  sr(x))  r+uk  +  J*  sr(s)APukds, 


wk+1  =  Twk  +  Buk,  k  =  0, 1 ,2,... 
w0eH 


where  T  e  £(H)  and  B  e  £(Rm,  H)  are  given  by  T  =  ST  (x)  and  B  =  (I  -  O'(x))  P  + 

JT  O'  (s)APds  respectively, 
o 

We  note  that  as  in  the  case  of  the  continuous-time  input  operator  B,  the  discrete- time  input 
operator  B  is  well  defined  and  does  not  depend  upon  a  particular  choice  for  P\  Indeed  if  Bj  and 
B2  are  the  input  operators  which  correspond  to  the  choices  P  and  r2+  then  for  ueRm  we  have 

-  b2)U  =  a  -  <r (x))(  r*  -  r2+)U  +  f  $r(s)A  (r/  -  r2+)uds. 


But  (P  -  P  )u  e  7l(T)  =  Dom(Ct)  and  therefore 


T  T 

J  $r(s)A(p  -  p  )uds  —  J  (T(s)Ct(p-p  )uds 

0  0 

T 

= «r (sxp  -  p  )uds  =  (^T  (x)  -  IXP-  -  rj )  u . 

0 

In  addition,  if  P  is  chosen  so  that  B,(P)  cz  71(A),  B  takes  on  the  particularly  simple  form  B  = 
(I  -  <T (x))P .  It  is  worth  noting  that  a  simple  calculation  reveals  that 


B  =  J  T  <T (s)  Bds 
o 


in  agreement  with  the  standard  technique  for  obtaining  the  discrete  or  sampled  time  formulation  of 
a  continuous  time  system  in  either  a  finite  dimensional  or  bounded  input  setting. 


It  is  our  intention  here  to  apply  the  approximation  theory  we  developed  earlier  in  (Gibson  and 
Rosen,  1986)  for  the  design  of  optimal  discrete-time  LQG  compensators  for  infinite  dimensional 
systems  with  bounded  input  and  output  operators.  We  therefore  require  the  additional  assumptions 
that  4)  T  =  tT (t)  e  Jo(V)  and  5)  ^(P)  c  V.  Although  not  all  boundary  control  systems  we 
might  formulate  would  satisfy  these  conditions,  there  are  many  interesting  and  important  systems 
which  do  (see,  for  example.  Section  4  below  and  Gibson  and  Rosen,  1987).  In  this  case,  the 
control  system  (2.1)  -  (2.4)  takes  the  form 


The  discrete-time  linear-quadratic  regulator  problem  for  the  boundary  control  system  (2.1)  - 
(2.3)  is: 

Find  u*  =  {u*k)k  =  0  e  -^(0,  °°;  Rm)  which  minimizes  the  quadratic  performance  index 
J(u)=  £<Qwk,wk>v  +  uIRuk 

fc= o 

where  Q  e  £(V)  is  self-adjoint  and  nonnegative,  R  is  a  symmetric  positive  definite  mxm  matrix 
and  the  state  w  =  (wk}~  Q  evolves  according  to  the  recurrence  (2.7),  (2.8). 

An  optimal  control  exists  for  each  initial  condition  Wq  if  and  only  if  the  operator  algebraic 
Riccati  equation 


(3.1) 


n  =  T*(n  -  riB(R  +  B*riB)IB*n)T  +  q  . 


'll] 


1 


t 


has  a  bounded  nonnegative,  self-adjoint  solution  n.  In  this  case,  the  optimal  control  has  the 

feedback  form  u^  =  -Fwk  where  F  =  (R  +  B*  riB)'1B*nT.  A  control  (sequence)  u  is  admissible 

for  the  initial  condition  wQ  if  the  corresponding  J(u)  is  finite.  If  there  exists  an  admissible  control 

for  each  initial  condition,  then  (3.1)  has  a  bounded  nonnegative,  self-adjoint  solution.  If  each 

admissible  control  for  each  initial  condition  drives  the  state  to  zero  asymptotically,  then  there  exists 

at  most  one  bounded  nonnegative,  self-adjoint  solution  to  (3.1).  The  optimal  trajectory  w*  = 

{w*v  }  evolves  according  to  w*  =  Skwn,  k  =  0,1,2,...,  where  the  closed  loop  state  transition 
k»0  k 

operator  S  e  JC(V)  is  S  =  T  -  BF.  If  Q  is  coercive,  then  S  has  spectral  radius  less  than  one  and  is 
uniformly  exponentially  stable.  From  the  finite  dimensionality  of  the  control  space  we  obtain 


uk  =  -<f,wk>v,  k  =  0,1,2,... 


where  f  =  (f1,f2,...,fm)T  6  x  Vis  called  the  optimal  functional  feedback  control  gain. 

i=i 

The  results  stated  here  for  the  optimal  linear-quadratic  regulator  problem  are  summarized  from 
(Gibson  and  Rosen,  1985). 

When  only  a  finite  dimensional  measurement  y  =  {yk}k  _  0  of  the  infinite  dimensional  state  w 
is  available  (recall  (2.9)),  a  state  estimator  or  observer  is  required.  For  a  given  input  sequence  u 
and  corresponding  output  sequence  y,  the  optimal  LQG  estimator  is 


(3.3)  wk+1  =  Twk  +  Buk  +  F  {yk  -  Cwk),  k  =  0,1,2,... 

(3.4)  w0  e  V 

A  A  A  ^  A  A  * 

where  the  optimal  estimator  or  observer  gain  F  e  £(Rp,V)  isF  =  TFIC  (R  +  dlC  )_1  with 

A* 

n  e  C(V)  the  minimal,  self-adjoint,  nonnegative  solution  (if  one  exists)  to  the  operator  algebraic 
Riccati  equation 


n  =  T(  n  -  fi  c*(R  +  c  fi  c Vc  n  )T*  +  q 


Since  F  e  (Rp,V),  it  has  the  representation 


-  T  P 

where  f  =  (f  t,  f  2,...,  fp) 1  £  x  V  is  called  the  optimal  functional  observer  gain. 

j=i 

A  A 

The  operator  Q  e  25(V)  is  self-adjoint,  nonnegative  and  the  pxp  matrix  R  is  symmetric, 
positive  definite. 

A  A 

In  a  stochastic  setting,  the  operator  Q  and  the  matrix  R  are,  respectively,  the  covariance 
operator  and  covariance  matrix  for  uncorrelated,  zero-mean,  stationary,  Gaussian  white  noise 

A 

processes  that  force  the  state  and  corrupt  the  measurement.  In  this  case,  if  Q  is  trace  class,  (3.3), 
(3.4)  is  the  infinite  dimensional  analog  of  the  discrete-time  Kalman-Bucy  filter.  In  a  strictly 

A  A 

deterministic  setting,  Q  and  R  are  assumed  to  be  determined  via  engineering  design  criteria 
such  as  stability  margins,  robustness  of  the  closed-loop  system,  etc. 

Replacing  operators  in  the  control  problem  with  the  adjoints  of  the  appropriate  operators  in  the 
estimator  problem  yields  the  usual  duality  between  the  LQG  optimal  control  and  estimator 
problems.  Hence  sufficient  conditions  for  existence  and  uniqueness  of  solutions  of  (3.5)  and  the 
closed-loop  estimator  stability  properties  are  analogous  to  the  results  for  the  control  problem.  In 
particular,  if  ek  =  w£  -  wk,  then  ek  =  £k  eQ,  k  =  0, 1,2, ...,  where  £  =  T  -  FC,  and  a  sufficient 
condition  for  £  to  be  uniformly  exponentially  stable  is  that  Q  be  coercive. 

The  optimal  LQG  compensator  consists  of  the  state  estimator  in  (3.3)  and  (3.4)  and  the  control 

law 

(3.6)  uk*  =  -F  wk*  >  k  =  0,1,2,... 

The  resulting  closed-loop  system  is  given  by 
V  k  =  /6k<Uf0,  k  =0,1,2,... 

where  Tlf  k  =  (wk,  w*k  )T  with  ( wk }  *k  =  0  the  state  trajectory  that  results  from  the  input  (3.6) 
and  A  e  C(VxV)  is 


S 

1 


•Sj 

•I* 

«! 

<2 


FC  T-BF-FC 


It  is  easy  to  show  that  the  spectrum  of  /4  is  given  by  a(/ 6)  =  a(S)  U  o(§),  so  that  the  stability  of 
the  closed-loop  plant-compensator  system  is  determined  by  the  stability  of  the  plant  with  full  state 
feedback  and  the  stability  of  the  estimator  error. 


For  each  N  =  1,2,...,  let  Vjsj  be  a  finite  dimensional  subspace  of  V  and  let  PN  be  a  bounded 
linear  mapping  from  V  onto  VN  (for  example,  the  orthogonal  projection  with  respect  to  either  the  V 

A 

or  H  inner  product).  Let  TN,  QN,  QN  e  X(VN),  BNe  £(Rm,VN)  and  CN  e  C(Vn,RP)  and  set 


Fn=  (R  +  B*nNBN)->B*nNTN 


FN=TNnNC*(R  +  CNhNCV 

A 

where  nN  and  IIN  are  the  minimal,  self-adjoint,  nonnegative  solutions  (assuming  that  they  exist)  to 
the  finite  dimensional  operator  algebraic  Riccati  equations 


nN  ~  TN(nN  *  nNBN(R  +  BHnNBN)-  B^nN)TN  +  qn 


-  tn(  nN  -  nN  (R  +  cNnN  cn)-'cn  nN)T*+  qn 


respectively.  The  approximating  optimal  compensator  is  given  by 


u*  =  -Fnw  ,k  =  0,1,2,.. 
H^X  N  H,X 


mm 


s 


> 


where  w  =  { w  }  “  is  determined  according  to  the  approximating  observer 
N  x-o 


w*  =  TN  w*  +  Bn  u*  +  FN{y*  -CNw*  ),  k  =  0,1,2... 
NX+1  N  N  n,x  n  -,m,x  n  nj 

w*  =  PN  w0  e  VN  . 

N,0  MON 


The  measurements  y  are  given  by  y  =  CwNJc  ,  k  =  0,1,2,...  where 

H  H  t  X 


wNJc+l=TwN.k-Bu*  k  •  k  =  0,1,2,... 


WN,0  =  w0 


The  resulting  closed-loop  system  is  given  by  <U/‘N  k  =  ^  *Uf  0,  k  =  0,1,2,...  where 

A  _ 

w 


Njc =  (wN.k-  w*  )T  and  ^NeZ  (VxVN)  is  given  by 

N  yX 


(3.9) 


/6n  .= 


-BF, 


N 


FC  T-BF-^C 
rN^  lN  DNrN  NT" 


Let  SN  =  Tn  -  BnFn  and  SN  =  TN  -  FNCN  and  assume  that  PN  ->  I  strongly  on  V  as 
N  -» <*».  Assume  further  that  TNPN  ->  T,  T*  PN  -»  T*,  QNPN  -4  Q  arid  QNPN  ->  Q 
strongly  on  V  and  that  BN  — *  B  and  CNPN  ~4  C  in  norm  as  N  4  ».  If  the  pairs  (TN,  BN) 
and  (T*  ,  C* )  are  uniformly  exponentially  stabilizable  and  the  pairs  (TN,  QN)  and  (T*  ,  QN)  are 
detectable  (see  Kwakemaak  and  Sivan,  1972)  then  there  exist  unique,  self-adjoint,  nonnegative 

A  A 

solutions  nNand  IIN  to  the  algebraic  Riccati  equations  (3.1)  and  (3.5).  If  nN  and  nN  are 

A  A 

bounded  from  above  uniformly  in  N,  then  nNPN  and  nNPN  converge  weakly  to  IT  and  n, 
respectively,  as  N->«  . 

A 

If,  in  addition,  SN  and  SN  are  uniformly  exponentially  stable,  uniformly  with  respect  to  N, 

A 

then  nNPN  and  I"INPN  converge  strongly.  Weak  convergence  of  nNPN  to  n  yields  strong 
convergence  of  FNPN  to  F  and  SNPN  to  S.  If  fINPN  converges  strongly  then  FNPN  -4  F  in 


norm.  Weak  convergence  of  nNPN  to  n  yields  weak  convergence  of  FN  to  F  and  SNPN  to 

A  A  A  A  A  A  A 

S.  When  nNPN  — »  FI  strongly,  then  FN  — >  F  in  norm  and  SNPN  — »  S  strongly  in  V  as 
N  ->  «».  Finally,  if  <PN  is  the  mapping  of  VxV  onto  Vx  ,'N  given  by  tPN(Wj,  w2)  =  (Wj, 

PNw2),  then  nNPN  — >  n  weakly  or  strongly  is  sufficient  to  conclude  that  — >  A  weakly 

/V  A 

or  strongly  depending  only  upon  whether  nNPN  — *  n  weakly  or  strongly  as  N  — » » .  Under 
appropriate  additional  hypotheses  on  the  spectral  properties  of  the  open-loop  system  and  on  the 
approximation  scheme,  it  is  possible  to  show  that  /6NiPN  converges  to  A  in  norm.  (We  have  been 
able  to  obtain  such  a  result  only  for  modal  approximations.)  Norm  convergence  of  the  closed-loop 
state  transition  operators  is  sufficient  to  conclude  that  uniform  exponential  stability  of  /(>  implies 
uniform  exponential  stability  of  for  all  N  sufficiently  large  (see  Gibson  and  Rosen  1986). 

In  practice,  the  finite  dimensional  approximating  subspaces  VN  are  often  constructed  using  any 
of  a  number  of  common  finite  element  bases,  e.g.  polynomial  and  hermite  spline  functions,  mode 
shapes,  orthogonal  polynomials,  etc.  For  the  discrete-time  boundary  control  systems  of  interest 
to  us  here,  the  approximations  to  T  and  B,  TN  and  BN,  are  obtained  by  approximating  the 
continuous  time  semigroup,  ( «T (t) :  t  £  0},  by  a  semigroup  of  bounded  linear  operators  on  VN, 
(TN(t):  t  £  0).  In  fact  it  is  the  infinitesimal  generator  Cl  of  the  semigroup  {ST(t):  t  >  0 } that  is 
approximated  by  a  bounded  linear  operator  ClN  on  VN  with  (TN(t):  t  >  0)  then  being  defined  by 
ayt)  =  exp  (CiNt),  t  >  0.  With  Tn=  Tn(t)  and  BN  =  (I  -  0yT))PNr+  +  r  0ys)PNAr+ds,  the 
required  convergence  can  usually  be  proved  using  the  Trotter-Kato  semigroup  approximation 

A  A 

result  (see  (Kato,  1966)  and  (Pazy,  1983) ).  The  approximations  to  Q,  Q  and  C,  QN  QN  and 

A  A 

CN ,  respectively,  typically  are  taken  to  be  QN  =  PNQ,  QN  =  PNQ  and  CN  =  CPN . 

Let{<pN}“N  denote  a  basis  for  VN  and  set  =  (cpN  ,cp^,...,  N*  ^n- 

Adopting  the  convention  that  [L]  denotes  the  matrix  representation  with  respect  to  the  basis 
(<p^)^  N  for  a  linear  operator  L  with  domain  and/or  range  in  VN,  we  find  that 
[Fn]  =  (R  +  [Bn]t  ©  N  [Bn])1[Bn]t  0  n[Tn]  and  [FN]  =  [TN]  ^[^^(R  +  [CNJ 

A  A 

0  n[Cn]t)'1  where  and  0  N  are  the  unique,  symmetric,  nonnegative  solutions  to  the 
nN  x  nN  matrix  algebraic  Riccati  equations 


n 


VM'il.'ikWil.ll.V.'ilMtiUiaU 


(3.10)  ©N  =  [Tn]t(©N-  ©^[BnKR  +  [Bn]T  0  n[Bn])1[Bn]t0n)[Tn] 

+  Mn(Qn] 


(3.1 1)  0N  =  [Tn](  ©N  -  Q^[Cu]t( R  +  [Cn]©n[Cn]t)-1[Cn]0n)[Tn]t 

+  [Qn](Mn)*1  . 

The  matrix  MN  is  the  nN  x  nN  Gram  matrix  «DN ,  (<DN)T>V  . 

If  w^k=(0N)TW^kw,tK  W*  ke  R^.then  G*  ^  -[FN)  ,  k  =  0,1,2,...  with 


A  A  A 


w;„r  tTN]W;jt+[BN]u-AhtFN](y^-[CN]W^).  k  =  0,!.2,. 


w*  =  (Mn)'1  <  <&N,  wn>v  . 

N,0  0  V 


111 

The  approximating  optimal  functional  feedback  control  gain,  f**  =  (f^,  f£...,  f^)T  e  x  VN 

j=i 

are  given  by  f**=  [FN](MN)'1<1>N  and  the  approximating  optimal  functional  observer  gain 

f*=(%  %...,  f$T  e  xPVN  by  f  =  [Fn]TOn.  IfnNPN  ->  n  weakly  (strongly) 

j=i 

A 

then  f1?  —>  fj,  i  =  1,2,...,  m  weakly  (strongly)  in  V.  If  nNPN  — »  n  weakly  (strongly)  then 
f1^  — >  fj  ,i  =  1,2,...,  p  weakly  (strongly)  in  V.  If  the  injection  V  <=  H  is  compact,  then  f^  -» f; , 

A  A  A 

i  =  1,2,...,  m  and  ft  -♦  fj.i  =  l,2,...,p  strongly  in  H  if  nNPN  and  FfNPN  converge  only 
weakly. 


We  consider  the  one-dimensional  heat  equation 


(4.1)  — (t,x)  =  a — —  (t,x),  0<x<l,  t>0. 


I 


K 


<• 


4. 


.  V  /  <"  «* 


i 

a 


»'t£/ 


I 


*jf  »tT  *»*  *■*  Cf  !.<  «,> 


»' ^UU  L.  L.  M -aai  *:.‘4V,  A  * «f. «.  A 


a 


where  a  >  0,  with  the  homogeneous  Dirichlet  boundary  condition  , 


(4.2)  w(t,0)  =  0,  t  >  0, 


and  either  the  Neumann  boundary  control 


(4.3)  — (t,l)  =  v(t),  t  >  0, 

dx 

or  the  Dirichlet  boundary  control 


(4.4)  w(t,  1)  =  v(t),  t  >  0, 

where  v  e  1^(0, »).  For  output  we  take  a  temperature  measurement 

(4.5)  y(t)  =  w(t,Q,  t  >  0, 

at  some  fixed  point  C,  e  (0,  1).  Initial  conditions  for  these  systems  have  the  form 

(4.6)  w(0,  x)  =  wQ  (x),  0  <  x  £  1 

where  w0  61^(0,  1). 

Although  the  two  control  systems  above  appear  to  be  similar,  they  are,  in  fact,  quite  different  and 
must  be  treated  separately.  We  begin  with  the  more  straight  forward  of  the  two,  Neumann  boundary 
control.  Let  H  =  Lj  (0,  1),  V  =  H^(0,  1)  =  {<p  e  H^O,  1) :  cp  (0)  =  0}  and 
W  =  H2  (0,  1)  fl  (0,  1).  With  H  endowed  with  the  usual  L2  inner  product,  V  with  the  inner 

product  <  (p,  V  >v  =  Jq  Dtp  D\p  and  W  with  the  inner  product  <  tp,  vp  >w  =  Jq  D)  Dd  \p,  we 

j=l 

have  the  continuous  and  dense  embeddings  WcVcHcV'c  W'.  Define  A  e  £( W,  H),  T  e  £(  W, 
R1)  by  A  tp  =  a  D2^,  Ttp  =  Dtp(l)  and  Ctp  =  tp(Q  respectively.  With  these  definitions  the  boundary 


.  .  i 

control  system  (4.1)  -  (4.3),  (4.5),  (4.6)  has  the  form  (2.1)  -  (2.4).  The  operator  Cl :  Dom(Ci)c 
H-»H  is  given  by  Cl(p=aD2<p  for  cp  e  {9  e  H2(0,  1):  9(0)  =  Dtp(l))  =  0}.  It  is  densely  defined, 
negative  definite,  self-adjoint  and  it  is  the  infintesimal  generator  of  a  uniformly  exponentially  stable  1 

analytic  semigroup  { <T (t):  t  >  0)  of  bounded,  self-adjoint  linear  operators  on  H.  Also,  { <T (t):  t  >  , 

0)  is  a  uniformly  exponentially  stable,  analytic  semigroup  of  bounded,  self-adjoint  operators  on  V  ] 


with  generator  Cj  given  by  0(p  =  Cup  for  tp  £  ftp  e  H3  (0, 1):  <p(0)  =  Dcp(l)  =  D2cp(0)  =  0}. 
Choosing  I'+  e  ^(R1,  W)  as  ( Ffu)(x)  =  xu  for  x  e  [0, 1],  we  have  1R1(  r+)  c  V,  <13»(  F4")  e  11(A) 
and  that  conditions  1)  -5)  given  in  Section  2  are  satisfied.  For  the  optimal  control  and  estimator 
problems,  we  take  Q  =  ql,  0  =  q  I,  R  =  r  and  ft.  =  ?  where  I  is  the  identity  on  V,  q,  q  >0  and  r,  r 
>  0.  The  uniform  exponential  stability  of  the  semigroup  { «T(t):  t  >  0}  on  V  implies  that  the  algebrai 
Riccati  equations  (3.1)  and  (3.5)  admit  unique  bounded,  nonnegative,  self-adjoint  solutions  n  and  fl 
respectively.  The  optimal  control  (3.2)  takes  the  form 


(4.7) 


DfDwfc,  k=  0,1,2,... 


where  the  optimal  functional  feedback  control  gain  f  and  the  optimal  functional  observer  gain 

A 

f  are  elements  in  Hl  (0,1). 

We  construct  an  approximation  scheme  using  a  linear  spline  based  Ritz-Galerkin  approach.  For 


each  N  =  1,2,... ,  (<pN}  N  q  denotes  the  usual  linear  spline  or  "hat"  functions  defined  on  the  interval 


[0,1]  with  respect  to  the  uniform  mesh  {0,  1/N,  2/N,...,l }.  We  discard  the  element  centered 
at  x  =  0,  tpN,  set  VN  =  span  {9^}^  ^  and  choose  PN  to  be  the  orthogonal  projection  of  V  onto  VN 
with  respect  to  the  V  inner  product.  Hence  VN  is  an  N  dimensional  subspace  of  V. 

For  <p  e  Dom  (0),  lCl<plH  >  al<plv  2:  al<plH  and  therefore  0  e  p(Ct)  and  CT*:  H  ->  Dom(Ct) 
satisfies  I CT^cpIy  ^  a''l<plH  for  9  e  H.  We  define  :  VN  — >  VN  as  the  inverse  of  the  operator 


=  PNCTl  restricted  to  VN.  The  operator -Ct^1  is  positive  definite  because 


<  Cf  1  9^  9n>v  =  -a*1l9Nl 2  for  9j^  e  VN ,  and  it  is  self-adjoint  since  <  0  "  1  9^  yN>v  = 

N  H  N 

<PnCT!  9n,  Yn>v  -  <<^’1  9n>  ¥n>V  =  a_1  <  9n*  Vn>h  •  Hence  the  operator  CiN  is  well  defined 
and  self-adjoint  For9N  eVNandyN  =CtN9N,  the  estimate 
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<  ^N(PN’(Pn>V  -  <  Vn,(-*n  HrN>V  -  "a  'l  Vn'h 

<  -alCrVN'v^  Vn’y  =  ‘al  aN‘ Wv 

=  _a  i  ^n'v 


implies  that  CiN  is  the  infinitesimal  generator  of  a  Gq  semigroup  { 3"N(t)  :  t  >  0}  of  bounded, 
self-adjoint  linear  operators  on  VN  satisfying  10"N(t)l  <,  e‘a  1  ,  t  >  0. 

It  can  be  shown  that  a«p  ,  \j/>v  =  <(-Ct)1^2<p,  (-Ci)^2y>H  .  It  then  follows  that  the  matrix 
representation  for  the  operator  CtN  with  respect  to  the  basis  {cp^} n  is  [ClN  ]  = 

-a<  cpN  <  (pN>  <pN>v  .  This  agrees  with  the  system  matrix  derived  by  a  standard 

Ritz-Galerkin  finite  element  approach.  Note  that  even  though  ClN  is  defined  to  be  the  inverse  of  the 
operator  PN(Ct)‘^  restricted  to  the  space  Vjq>  computing  its  matrix  representation  does  not  require 
either  Ct 1/2  or  CT1  explicitly.  In  general,  the  same  approach  can  be  used  to  obtain  an  operator 
representation  for  the  Ritz  Galerkin  approximation  to  any  self-adjoint  coercive  operator. 

Let  SJN  denote  the  interpolation  operator  from  V  onto  VN  defined  by  (dN  <p)(j/N)  =  <p(j/N),  j  =  1, 
2,...,  N  .  Then  for  tpeW,  elementary  approximation  properties  of  linear  interpolator  spline 
functions  (see  (Schultz,  1971))  imply 


l(PN  -  I)<plv  <;  l(dN  -  I)cply  ^  —  !D29(h 

N7t 

and  therefore,  since  W  is  dense  in  V,  that  PN  -»  I  strongly  on  V  as  N— >  °°  .  Also,  it  follows  that 
Cl^  =  PNCr^  Cl'1  strongly  on  V  as  N->  oo  .  If  we  define  TN  =  <Tn(t),  then  the  Trotter-Kato 
approximation  theorem  yields  that  TN  Pn-^T  strongly  on  V  as  N— >  °°  and,  since  T  =  T  =  «T(t) 
and  =  TN  =  tT j^x)  that  T^,Pn  T*  strongly  on  V  as  N->  «>. 

Since  fflOT*  )cvN  (recall  that  (r^uXx)  =  ux,  0  ^  x  <  1),  we  define  the  approximating  input 

A 

operators  BN  by  BN  =  (I  -  tTN(x))P1'  and  set  QN  =  ql,  QN  =  q  I  and  CN  =  C.  The  strong  convergence 
of  PN  to  the  identity  and  tnpn  to  T  together  with  the  finite  dimensionality  of  the  domain  of  B  and  the 

A  A 

range  of  C  are  sufficient  to  conclude  that  QNPN  -» Q,  QNPN  — »  Q  strongly  on  V  and  that  BN  — »  B 
and  CNPN  -4  C  in  norm  as  N  -»  «>. 

The  uniform  exponential  stability  of  the  semigroups  {  iTN(t)  :  t  >  0}  implies 


, 


'd 

VJ 

'Ijj 

!t1 


TO 


.lil 


I 


(4.8)  ITk  iv  =  l(T*  )kiv  <  r*  k  =  0, 1 ,2,... 

N  v  N  v 


with  r  =  e  at  <  1 .  Consequently  the  pairs  (TN.BN)  and  (TN,  C* )  are  uniformly  exponentially 

A 

stabilizable  and  the  pairs  (TN,QN)  and  (TN,  QN)  are  detectable.  It  follows  that  there  exist  unique 

A 

self-adjoint,  nonnegative  solutions  nN  and  nN  to  the  finite  dimensional  algebraic  Riccati 
equations  (3.7)  and  (3.8)  respectively.  The  uniform  exponential  bound  (4.8)  with  r  <  1  imples  that 

A 

the  zero  control  yields  a  uniform  upper  bound  for  nN  and  nN  and  therefore  the  uniform 

A  A 

exponential  stability  of  SN  =  TN  -  BNFN  and  SN  “  TN  "  FNCN  •  We  conclude  that  nNPN  and 

A  A  Aw 

nNPN  converge  strongly  in  V  to  nN  and  nN  ,  respectively  ,and  that  FNPN  and  Fn  converge  to  F 

A 

and  F  in  norm  as  N  — »  The  approximating  optimal  functional  feedback  control  and  observer 

A  A 

gains,  fN  and  fN,  converge  respectively  to  f  and  f  in  the  H1  norm  as  N  — >  «>. 

In  implementing  the  approximation  scheme  just  outlined  above,  eigenvector  decomposition  of 
the  associated  Hamiltonian  matrix  was  used  to  solve  the  matrix  algebraic  Riccati  equations  (3.10)  and 
(3.1 1)  (see  Pappas,  et.  al.,  1980).  The  required  matrix  exponentials  also  were  computed  using 
eigenvalue/eigenvector  decomposition.  All  calculations  were  carried  out  via  Fortran  codes  on  an 

_  A  __ 

IBM  PC  AT.  We  set  a  =  V.l,  q  =  q  =  r  =  r  =  1.0,  5,  =  V2/2  and  t  =  .01  and  obtained  the 

A  A 

functional  gains  plotted  in  Figs.  4.1  and  4.2  .  We  plot  fN  and  fN  as  well  as  DfN  and  DfN  to  exhibit 
the  H1  convergence.  We  note  that  Df  (or  DfN)  appears  as  the  feedback  kernel  in  the  optimal  control 
law  (4.7). 

We  also  simulated  the  operation  of  the  closed-loop  system  with  an  approximating  compensator. 
Using  a  20  mode  model  for  the  infinite  dimensional  system  and  N  =  12,  we  computed  the  closed- 
loop  spectrum  of  the  approximating  compensator  (i.e.  the  eigenvalues  of  the  operator  /6N  given  by 
(3.9)  with  N  =  12).  These  eigenvalues  along  with  the  first  20  open-loop  eigenvalues  (i.e.  the  first 
20  eigenvalues  of  the  operator  T  =  <T(x))  and  the  approximating  closed- loop  control  and  observer 
eigenvalues  are  tabulated  in  Table  4.1  below.  Table  4.1  reveals  that  the  last  seven  open-loop 
eigenvalues  remain  essentially  unchanged  in  the  closed-loop  system-i.e.  these  modes  are  neither 
controlled  nor  observed  by  the  finite  dimensional  compensator.  Also,  as  one  would  expect,  a(^6N) 

A 

consists  essentially  of  the  union  of  o(SN) ,  <J(SN)  and  the  eigenvalues  corresponding  to  the 


uncontrolled/unobserved  modes  of  the  open-loop  system. 

It  is  worth  noting  that  the  scheme  we  have  outlined  above  for  the  Neumann  boundary  control 
problem  is  the  same  scheme  that  one  would  ordinarily  use  if  the  problem  were  formulated  in  the  space 
H  -  i.e.  if  the  output  operator  C  was  bounded  on  11^(0, 1)  (see  Gibson  and  Rosen,  1986).  This  is 
possible  primarily  because  the  space  V  =  H^(0,1)  is  the  natural  energy  space  for  the  underlying 
homogeneous  or  open-loop  system.  Consequently,  the  inherent  self-adjointness  and  coercivity  in  the 
problem  is  preserved  when  it  is  formulated  in  the  stronger  space.  In  the  case  of  Dirichlet  boundary 
control,  the  situation  is  quite  different. 

For  the  Dirichlet  boundary  control  system  (4.1),  (4.2),  (4.4)  -  (4.6),  we  choose  the  spaces  H,  V 

and  W  and  their  corresponding  inner  products  to  be  the  same  as  they  were  in  the  Neumann  case. 

The  operators  A  £  2S(W,H)  and  C  e  ^(V.R1)  also  remain  unchanged,  however  now  we  have 

T  e  X>(W,Rl)  given  by  Tcp  =  cp(l).  It  then  follows  that  the  operator  Ct  :  Dom(Ct)  c  H  — >  H  is  given 

by  Crip  =  aD2<p  for  tp  e  H2(0,1)  D  H1  (0,1).  It  is  well  known  that  Ct  is  densely  defined,  negative 

o 

definite  and  self-adjoint  and  that  it  is  the  infinitesimal  generator  of  the  uniformly  exponentially  stable 

analytic  semigroup  { S' (t) :  t  >  0}  of  bounded,  self-adjoint  linear  operators  on  H.  However  this 

time  the  operators  S' (t)  for  t  >  0  are  neither  self-adjoint  nor  a  semigroup  on  V.  Indeed,  since 

fll(<T(t))c  H1  (0,1)  for  all  t  >  0  and  H1  (0,1)  is  a  closed  proper  subspace  of  H1  (0,1),  S(t)  is  not 
o  o  L 

strongly  continuous  in  the  V-norm  at  t  =  0.  (The  fact  that  our  general  framework  requires 

TP'  =  1  and  (R.fr4')  c  V  precludes  our  choosing  V  to  be  H1  (0,1).)  On  the  other  hand, 

o 

(S(t) :  t>  0}  an  analytic  semigroup  implies  (see  Pazy,  1983)  that  there  exists  a  constant  p.>  0  foi 
which  lets' (t)lH  <  |0.f  *  for  t  >  0  .  Consequently,  if  we  define  T  =  S' (x),  then  it  follows  that 
T  e  £(V)  and  moreover,  that 


ITkcpl2 

v 


=  -  a'^<CtS' (kx)cp.S' (kx)cp>  <  a'1 !  QS(kx)tplH  (S' (kx)cp|H 


-akt  -akt 

ue  ,  ,2  .  pie  2 
£  - Itpl  <  - lcplv 

akx  akx 


VV  *  *  • 


for  k  =  1,2,...  and  <p  e  V.  We  have  therefore 


(4.9)  ITklv  =  l(T*)klv  <;  Mrk,  k  =  0,1,2,... 


where  M  >  0  and  r  <  1. 

We  again  choose  I'4'  e  jGCRSW)  as  (  T'hjXx)  =  xu  for  x  e  [0,1].  Then  R(P1')  c  Tl(A)  and  we 
have  reformulated  the  boundary  control  system  (4.1),  (4.2),  (4.4)  -  (4.6)  in  the  general  form  of  (2.1) 
-  (2.4)  and  conditions  1)  -  5)  are  satisfied. 

We  formulate  the  optimal  control  problem  with  the  performance  index 

J(u)=  £  q<wk,wk>H  +  ru^ 
k=0 


where  q  >  0  and  r  >  0.  That  is,  we  take  Q  to  be  the  bounded,  self-adjoint  nonnegative  operator  on 

1  X  1 

Hl(0,1)  given  by  (Qcp)(x)  =  q  jQ  J  tp( z)dzdy  and  R  to  be  r.  For  the  estimator  problem  we  set 

A  A  A  A 

Q  =  q  I  and  R  =  r  with  q  >  0  and  ?  >  0. 

The  uniform  exponential  bound  (4.9  )  implies  the  existence  of  unique,  nonnegative,  self-adjoint 

A 

solutions  FI  and  fl  to  the  algebraic  Riccati  equations  (3.1)  and  (3.5).  The  optimal  control  is  again  of 


the  form  (4.7)  with  the  optimal  functional  gains  f  and  f  in  HL . 

The  fact  that  { (T(t) :  t  >  0}  is  not  a  semigroup  on  V  precludes  the  use  of  a  semigroup  -  theoretic 
approach  to  approximation.  We  therefore  employ  modal  subspaces  and  approximate  the  open-loop 
state  transition  operator  T  directly  as  a  bounded  linear  operator  on  V. 

For  each  N  =  1,2,...  let  VN  =  span  {<pj }  ^  where  for  x  e  [0,1],  cp0(x)  =  x  and  tpj(x)  =  sinjnx, 

j  =  1,2,. ...N.  Let  pN  denote  the  orthogonal  projection  of  H  =  L2(0,1)  onto  span  {cpj }  N 


and  let  PN  denote  the  orthogonal  projection  of  V  onto  VN.  Using  the  fact  that  V  =  FU(0,1)  ©  (pg>  ^ 


is  not  difficult  to  see  that  PNcp  =  <p(l)<p0  +  pN(cp  -  cp(I)cp0)  for  cp  e  V  and  hence,  via  elementary 
properties  of  Fourier  series  (see  Tolstov,  1962),  that  I(PN  -  r)(p|y  =  l(pN  -  I)(cp  -  cp(l)<p0)|v  ->  0 
as  N  — »  oo  for  each  <p  e  V. 

We  define  TNe£(VN)  by  Tn  =  PnT.  Then,  since  fi.(T)  =  R(ST(t))  c  H‘(0,1), 
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$ 

I 
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for  VN  =  X  Vn  ‘Pj  ■ 6  VN  we  have 
j=o 

^  2(  lV  0  ^Tl^T 

tn  Vn  =  PNTyN  =  PN^r (t)VN  =PN^ (t)v  =  T (t)pNvN  =  2^  {  VN  +  VN}  e' aJ  *  \  . 

j=i  J71 

It  follows  that  =  PNT*,  iTk  ly  =  l(T^)klv  <  MrK  k  =  0,1,2,...  with  M  >  0  and  r  <  1  independent 
of  N,  and  that 


l(TNPN  -  T)cplv  <  l(PNTPN  -  PNT)cptv  +  l(PN  -  I)Tcplv 
<  Mrl(PN  -  I)(plv  +  l(PN  -  I)Tcp|v  ->  0 


as  N  ->  <»  for  tp  e  V.  Similarly,  T'  ‘  PN  ->  T*  strongly  on  V  as  N  ->  <». 

A 

The  approximating  input,  output,  and  state  penalization  operators  BN,  CN,  QN  and  QN  take  the 


Bnu  =  (I  -  TN)r*u  =  <p0u  +  £  2^  e  9ju  , 


j=l  J71 

A 

CN  ~C  ,Qn  =  qPNQ  and  QN  =  q  I.  Reasoning  as  we  did  in  the  Neumann  case,  the  approximating 
algebraic  Riccati  equations  (3.7)  and  (3.8)  admit  unique,  nonnegative,  self-adjoint  solutions  nN  and 

A  A  A  /N  A 

nN  respectively,  nNPN  — >  FI  and  nNPN  II  strongly  on  V  and  FNPN  -»  F  and  FN  F  in  norm 

as  N  — >  «> .  The  approximating  functional  feedback  control  and  observer  gains  fN  and  fN  converge 

A 

to  f  and  f  respectively,  strongly  in  H1  as  N  — »  «> . 

With  a  =  1.0,  q  =  q  =  r  =  1.0,  ?  =  5.0,  ^  =  V2/2  and  %  =  .01  and  the  scheme  oudined  above  we 
obtained  the  approximating  optimal  functional  feedback  control  and  observer  gains  plotted  in  Figs. 
4.3  and  4.4  below.  The  first  12  open-loop  and  the  approximating  closed-loop  control  and  observer 
eigenvalues  for  N  =  12  are  tabulated  in  Table  4.2. 

Table  4.2  reveals  an  interlacing  of  the  closed-loop  control  and  open-loop  eigenvalues.  That  is, 
the  closed-loop  control  eigenvalues  (i.e.  the  elements  in  the  spectrum  of  S)  are  alternately  more  and 
less  stable  than  the  corresponding  open-loop  eigenvalues.  We  also  have  observed  this  phenomenon 
in  other  numerical  studies  we  are  carrying  out  involving  LQG  boundary  control  for  flexible 
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•  structures.  In  additon,  in  the  Dirichlet  boundary  control  system  discussed  above,  if  Q  is  chosen  as 
the  identity  operator  on  V  =  H[(0,1),  virtually  all  of  the  closed-loop  control  eigenvalues  are  less 
stable  than  the  corresponding  open-loop  eigenvalues.  It  is  clear  that  this  non-standard  behavior 
results  from  the  presence  of  the  one  dimensional  subspace  represented  by  <H(rf).  Indeed,  the 
behavior  of  the  closed-loop  spectrum  in  the  case  of  Neumann  boundary  control  is  as  would  be 
expected.  We  feel  that  what  we  are  seeing  can  most  likely  be  explained  via  infinite  dimensional 
analogs  of  existing  results  relating  the  asymptotic  properties  of  the  closed-loop  spectrum  of  a  linear 
regulator  and  the  zeros  of  the  corresponding  open-loop  transfer  function  (see  Kwakemaak  and  Sivan, 
1972  and  Harvey  and  Stein,  1978).  However,  as  of  yet,  we  have  been  unable  to  establish  this 
conjecture  satisfactorily  and  we  consider  it  to  be  beyond  the  scope  of  this  paper,  which  is  primarily 
concerned  with  approximation.  We  leave  it  as  an  interesting  open  question. 


We  have  developed  a  framework  for  the  finite  dimensional  approximation  of  optimal  discrete-time 
LQG  compensators  for  distributed  parameter  systems  with  boundary  input  and  unbounded 
measurement.  Our  theory  applies  to  the  class  of  boundary  control  problems  which  can  be  formulated 
in  a  state  space  in  which  both  the  discrete-time  input  and  output  operators  are  continuous.  We  have 
used  a  functional  analytic  treatment  to  develop  a  convergence  theory  and  have  demonstrated  the 
feasibility  of  our  approach  via  examples  involving  either  the  Neumann  or  Dirichlet  boundary  control  of 
a  one  dimensional  heat  equation  with  point  measurement  of  temperature.  We  have  shown  that  while 
both  problems  outwardly  appear  to  be  quite  similar,  they  in  fact  require  very  different  approaches  to 
approximation.  Also  in  the  Dirichlet  case  the  observed  behavior  of  the  resulting  closed-loop  spectrum 
is,  in  some  ways  unexpected  and  its  explanation  remains  open. 

Finally,  we  have  been  looking  at  the  application  of  our  schemes  to  LQG  problems  for  flexible 
structures  with  boundary  inputs  and  unbounded  measurement  and  systems  with  control  and/or 
observations  delays.  We  have  been  considering  vibration  suppression  for  cantilevered  beams  via 
shear  or  moment  inputs  at  the  free  end  and  pointwise  observation  of  strain  or  acceleration.  These 
studies  are  currently  underway  with  the  results  to  be  reported  elsewhere. 


$ 

% 

i 

l 


Jm*  A* k* .h* A’ th<A'  Jt*  J|*  J»*  |h*  t*’>f  ll^  W IMCy  ll‘  W  h4*i<t * **■  •* »*«>.  *»>  i^nu  tWIA  -kL. ~ il  -« 


Open-Loop 


oC*12) 


°(s 12) 


o(S12) 


.9975 

.9968 

.9968 

.9780 

.9780 

.9778  | 

.9769 

.9768 

.9402 

.9408  ! 

.9387 

.9371 

.9371 

.8861 

.8872 

.8778 

.8775 

.8798 

.8188 

.8194 

i 

.7982 

.7985 

.7998 

.7419 

.7414 

.7026 

.7019 

.7030 

.6590 

.6573 

.5960 

.5921 

.5946 

.5740 

.5718 

.4891 

.4769 

.4804 

.4901 

.4875 

.4433 

.4412 

.4104 

.4041 

.3675 

.3675 

.3682 

.3368 

.3341 

.2772 

.2763 

.2768 

.2711 

.2705 

.2145 

.2129 

.2133 

.2139 

.2134 

.1811 

.1811 

.1816 

.1655 

.1663 

.1255 

.1260 

.0934 

.0933 

.0681 

.0677 

.0482 

.0483 

.0341 

.0340 

.0235 

.0236 

‘iWi1 


)>SviVjvivV» 


'p’». 


s< 


A 

Open-Loop 

o(S12) 

°(s 12> 

.90601806 

.67382545 

.41136911 

.20615299 

.08480497 

.02863695 

.00793790 

.00180617 

.00033753 

.00005172 

.00000651 

.00000067 


.90569591 

.78573771 

.68243047 

.57981918 

.40961171 

.40082268 

.20758391 

.20936323 

.08447005 

.08636884 

.02873534 

.02892353 

.00791793 

.00792193 

.00180978 

.00178763 

.00033682 

.00033414 

.00005179 

.00005162 

.00000650 

.00000654 

.00000067 

.00000068 

.00000000 

.00000000 

MWtl 


. 

Il 


1$$: 


Curtain,  R.  F.  (1984).  Finite  dimensional  compensators  for  parabolic  distributed  systems  with 
unbounded  control  and  observation.  SIAM  J.  Control  and  Opt.,  22,  pp.  255-276. 

Curtain,  R.  F.  and  D.  Salamon  (1986).  Finite  dimensional  compensators  for  infinite  dimensional 
systems  with  unbounded  input  operators.  SIAM  J.  Control  and  Opt.,  24,  pp.  797-816. 

Gibson,  J.  S.  and  I.  G.  Rosen  (1985).  Numerical  approximation  for  the  infinite  dimensional 
discrete-time  optimal  linear  quadratic  regulator  problem.  ICSAE  Report  No.  86-15,  Institute  for 
Computer  Applications  in  Science  and  Engineering,  NASA  Langley  Research  Center,  Hampton,  VA, 
and  SIAM  J.  Control  and  Opt.,  to  appear. 

Gibson,  J.  S.  and  I.  G.  Rosen  (1986).  Computational  methods  for  optimal  linear  quadratic 
compensators  for  infinite  dimensional  discrete-time  systems.  ICASE  Report  No.  86-72,  Institute  for 
Computer  Applications  in  Science  and  Engineering,  NASA  Langley  Research  Center,  Hampton,  VA, 
and  Proceedings  of  Conf.  on  Control  and  Identification  of  Distributed  Parameter  Systems,  Vorau, 
Austria,  July,  1986,  to  appear  as  a  volume  in  the  Springer- Verlag  Lecture  Notes  in  Control  and 
Information  Sciences. 

Gibson,  J.S.  and  I.G.  Rosen  (1987).  Approximation  in  discrete-time  boundary  control  of  flexible 
structures.  Proceedings  of  The  26th  IEEE  Conference  on  Decision  and  Control,  December  9-11, 
1987,  Los  Angeles,  California. 

Harvey,  C.A.  and  G.  Stein  (1978).  Quadratic  weights  for  asymptotic  regulator  properties. 

IEEE  Trans.  Auto.  Control,  AC-23,  pp.  378-387. 


Hille,  E.  and  R.S.  Phillips  (1957).  Functional  Analysis  anc 
Society,  Providence,  RI. 

Kato,  T.  (1966).  Perturbation  Theory  for  Linear  Opera  ton 


ri -Groups.  American  Mathematical 


Springer-Verlag,  New  York. 


Kwakemaak,  H.  and  R.  Sivan  (1972).  Linear  Optimal  Control  Systems.  Wiley-Interscience,  New 
York. 

Lasiecka,  I.  and  R.  Triggiani  (1983a).  Dirichlet  boundary  control  problem  for  parabolic  equations 
with  quadratic  cost:  analyticity  and  Riccati  feedback  synthesis.  SIAM  J.  Control  and  Opt ,  21  dd 
41-67.  F  FF 

Lasiecka,  I.  and  R.  Triggiani  (1983b).  An  L,  theory  for  the  quadratic  optimal  cost  problem  of 
hyperbolic  equations  with  control  in  the  Dirichlet  boundary  conditions.  In  F.  Kappel, 

K.  Kunisch  and  W.  Schappacher  (Eds.),  Control  Theory  for  Distributed  Parameter  Systems  and 
Applications,  Lecture  Notes  in  Control  and  Information  Sciences  54,  Springer-Verlag  Berlin  pd 
138-152.  e  FF' 

Lasiecka,  I.  and  R.  Triggiani  (1986).  Riccati  equations  for  hyperbolic  partial  differential  equations 
with  L2(0,  T;  Lj  (D)  -  Dirichlet  boundary  terms.  SIAM  J.  Control  and  Opt.,  24,  pp.  884-925. 

Lasiecka,  I.  and  R.  Triggiani  (1987a).  The  regulator  problem  for  parabolic  equations  with  Dirichlet 
boundary  control  Part  I:  Riccati’s  synthesis  and  regularity  of  optimal  solution.  J.  Appl.  Math,  and 
Opt.,  to  appear. 

Lasiecka,  I.  and  R.  Triggiani  (1987b).  The  regulator  problem  for  parabolic  equations  with  Dirichlet 
boundary  control  Part  II:  Galerkin  approximation.  J.  Appl.  Math,  and  Opt.,  to  appear. 

Pappas,  T.,  A.  J.  Laub  and  N.  R.  Sandell,  Jr.  (1980).  On  the  numerical  solution  of  the  discrete-time 
algebraic  Riccati  equation.  IEEE  Trans.  Auto.  Control,  AC-25,  pp.  631-641. 


APPROXIMATION  IN  DISCRETE-TIME  BOUNDARY  CONTROL  OF  FLEXIBLE  STRUCTURES 


1 


I * 
#| 


a 


JS.  Gibson 

Department  of  Mechanical,  Aerospace  and  Nuclear  Engineering 
University  of  California 
Los  Angeles,  CA  90024-1597 

I.G.  Rosen 

Department  of  Mathematics 
University  of  Southern  California 
Los  Angeles.  CA  90089 


ABSTRACT 

This  paper  treats  discrete-time  LQG  optimal 
control  of  flexible  structures  with  boundary 
control  and  Mhat  normally  are  unbounded  measurement 
operators.  The  application  of  recently  developed 
approximation  theory  for  Infinite  dimensional 
discrete-time  LQG  problems  to  the  problem  here  is 
discussed,  and  a  numerical  example  is  presented. 

1 .  INTRODUCTION 

Recently,  approximation  theory  has  been  deve¬ 
loped  for  the  linear-quadra^ ic-Guasslan  (LQG) 
control  problem  for  infinite  dimensional  systems 
with  bounded  input  and  output  operators,  and  the 
application  for  optimal  digital  control  of  flexible 
structures,  diffusion  processes  and  time-delay 
systems  has  been  demonstrated  in  numerical  examples 
(see  [3.4.5]).  In  the  corresponding  continuous- 
time  control  problem,  the  requirement  that  the 
input  operator  be  a  bounded  operator  from  the 
finite  dimensional  control  space  into  the  Infinite 
dimensional  state  apace  almost  always  eliminates 
boundary  control  of  distributed  systems.  However, 
in  discrete-time  control  with  piecewise  constant 
inputs,  the  input  operator  representing  a  boundary 
control  usually  la  bounded.  This  makes  the  theory 
developed  in  [3]  and  [4]  applicable  to  boundary 
control  of  distributed  models  of  flexible  struc¬ 
tures.  An  important  feature  of  the  approach  in 
this  paper  is  that  we  accommodate  both  boundary 
control  and  a  normally  unbounded  measurement  opera¬ 
tor  by  working  the  problem  in  a  Hilbert  space  (the 
space  V  in  the  following  sections)  with  stronger 
topology  than  the  natural  energy  space  (the  space  X 
in  the  following  sections).  This  trick  appears  to 
simplify  many,  though  not  all,  problems  whose 
continuous-time  analogs  are  much  more  difficult. 

In  this  paper  we  sketch  the  theory  for  trans¬ 
forming  an  abstract  boundary  control  system  with 
unbounded  measurement  operator  into  a  discrete-time 
control  system  with  bounded  input  and  output  opera¬ 
tors.  Also,  we  summarize  the  discrete-time  LQG 
approximation  theory  in  [3]  and  14],  and  we  give 
numerical  results  for  control  of  a  flexible  bean. 

2 .  THE  ABSTRACT  BOUNDARY  CONTROL  PROBLEM 
AND  THE  DISCRETE-TIME  FORMULATION 

We  consider  boundary  control  problems  for 
flexible  structures  whose  equations  of  motion,  in 
first-order  form,  can  be  written  abstractly  as 
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x ( t )  =  Ax ( t ) ,  t  >  0,  (2.1) 

Tx(t)  -  u ( t ) ,  t  >  0.  (2.2) 

X (0 )  -  xQ,  (2.3) 

y ( t )  -  C  x(t),  t  >  0,  (2.4) 

where  A  e  L(W,X),  T  e  L(W.r“),  C  c  L(V.RP)  and 
xQ  €  X  with  X,  W  and  V  Hilbert  spaces.  We  assume 
that  W  and  V  are  densely  and  continuously  embedded 
in  X;  that  the  operator  T  Is  surjective  and  its  null 
space.  N(T),  is  dense  in  X;  and  that  the  operator 
A : Don ( A )  c  X  -  X.  defined  to  be  the  restriction  of 
A  to  NIT),  generates  an  CQ-semigroup  (T(t):  t  >  0) 
of  bounded  linear  operators  on  X. 


We  treat  the  unbound  input  by  rewriting  the 
boundary  control  system  (2.1) -(2. 3)  as  an  equiva¬ 
lent  evolution  system  in  a  space  Z,  larger  than  X, 
and  then  considering  weak,  or  mild,  solutions  (see 
[2]  for  this  approach).  To  construct  Z,  we  begin 
with  the  space  Dom(A*)  endowed  with  the  graph 
Hilbert  space  norm  associated  with  A*.  We  then 
define  the  dual  of  Dom(A*)  to  be  Z.  It  can  be 
shown  that  X  is  densely  and  continuously  embedded 
in  Z  and  that  (T(t):  t  <  0)  admits  a  unique  exten¬ 
sion  to  a  Cg-semlgroup  on  Z.  The  generator  of  this 
extended  semigroup  is  an  extension  A  of  A.  with 
A:X  c  Z  -  Z  given  by  ( Ad> ) (|»  -  <$,A*i|i>  for  <f  e  X  and 
if  €  Dorn  (A* )  . 

We  recall  that  T  was  assumed  to  be  onto  and 
let  T*  e  L(R*,W)  denote  any  right  inverse  of  T. 

Then  for  u  e  R*.  we  define 


§u  -  AT*u  -  ArV 


It  is  not  difficult  to  show  that  B  is  a  well  defined 
element  of  L(R.Z),  which  does  not  depend  on  the 
partiuclar  choice  of  T+.  For  u  e  L2(o,tf:  r")  and 
x0  e  X.  it  can  be  shown  that  the  function 


x(t)  -  T(t)  xQ  ♦  f  T(t-s)Bu(s)ds 
0J 


Is  in  C([0.tf];X)  n  H  ( 0 ,  t  ;  Z ) ,  that  x(0)  »  xQ  and 

X ( t )  =■  A  x(t)  +  Bu(  t )  (2.7) 

In  Z  for  alaost  every  t  e  [0,t^].  The  x-valued 
function  x  given  by  (2.6)  is  called  a  weak  solution 
to  the  Initial  value  problem  with  boundary  input, 
(2. 1 ) - (2 . 3)  . 


We  derive  the  discrete-tine  formulation  of 
(2.1)— (2.3)  by  letting  T  denote  the  length  of  the 
sampling  interval  and  considering  piecewise 
constant  control  Inputs  of  the  form  u(t)  «  u(kr)  , 

kT  <  t  <  (k+1 )T.  k  *  0,  1 ,  2 .  From  (2.6) 

then,  we  obtain 


x((k+l)T)  *  T  x(kr)  ♦  B  u(kt) 


T  -  T(kr)  e  L(X) , 


-T 

b  =  ( i  -  t(t) >r"  *  T(t)ar+dt  c  l(r”,x>. 


Note  that  if  r  can  be  chosen  so  that  R(T  )  c  N( A) 
then  the  discrete-tine  input  operator  B  takes  the 
simple  form  B  ■  (I  -  T)f+. 


To  simplify  notation,  we  will  assume  henceforth 
that  the  tine  scale  has  been  normalized  to  make 
T  =  1 .  The  discrete-t ime  control  system  Is  then 


x(k*l)  -  T  x(k)  *  B  u(k) 


k  ■=  0 ,  1,  2 . 


y(k)  =  C  x(k) . 


If  the  Initial  state  x(0)  =  xQ  e  X  as  we  originally 
assumed,  then  each  x(k)  in  (2.10)  will  be  in  X,  so 
that  the  discrete-time  control  system  does  not 
involve  the  space  Z. 


From  here  on.  we  will  assume  that  x ( 0 )  is  In 
the  stronger  space  V.  The  reason  is  that  we  will 
apply  the  approximation  theory  developed  in  [3]  and 
(4]  for  optimal  discrete-time  LQG  compensator 
synthesis  for  systems  with  bounded  Input  and 
bounded  output  operators  to  boundary  control 
systems  with  measurement  operators  that  are 
unbounded  on  the  a pace  X.  While  the  discrete-time 
input  operator  B  is  an  element  in  1. ( R . X ) .  in 
general  the  output  operator  C  Is  not  bounded  from  X 


to  R  .  Consequently,  we  must  consider  the 
discrete-time  control  system  ( 2 . 10) - (2 . 11 )  in  the 
space  V,  where  both  the  input  and  output  operators 
are  bounded.  In  addition  to  x(0)  c  V,  therefore, 
we  require  that  T  -  T(r)  £  L(V)  and  R(T)  c  v,  and 
we  observe  that  these  conditions  imply  B  e  L(R,V) 
and  x(k)  c  V,  k  -  0,  1,  2,  ...  .  While  not  all 
boundary  control  problems  for  flexible  structures 
satisfy  these  hypotheses,  a  large  class  of  impor¬ 
tant  problems  do.  This  should  be  evident  from  the 
example  in  Section  4. 


3.  THE  0 PTiWIAL  DISCRETE-TIME  LQG 
COMPENSATOR  AND  APPROXIMATION 
We  present  here  a  brief  summary  of  the  linear- 
quadratic  theory  and  associated  approximation 
results  developed  in  (3],  [4]  and  [5]  for  Infinite 
dimensional  discrete-time  control  systems  having 
the  form  (2 . 10 )- (2 . 11 )  ,  with  T  e  L(V),  B  £  L(R*.  V) 

p 

C  e  L(V ,  R  (and  x(0)  e  V.  The  optimal  compensator 
is  an  optimal  feed  back  control  law  together  with 
an  optimal  state  estimator,  or  observer.  We  treat 
the  control  problem  first. 


The  linear-quadratic  optimal  control  problem  is 
to  find  a  control  sequence  u ( • )  to  minimize  the 
performance  index 


J(u)  =  Y_  <  Q  x,k)-  x,k)  >  *-u<k>T  R  u(k)  (3.1) 


where  x(0)  £  V,  x(k).  k  *  1,  2 . is  given  by 

(2.10),  Q  £  L(V)  is  nonnegative  and  self  adjoint, 
and  R  is  an  m  x  m  positive  definite  symmetric 
matrix.  For  this  problem  the  infinite  dimensional 
theory  closely  parallels  the  finite  dimensional 
case.  A  solution  to  the  optimal  control  problem 
exists  for  each  x ( 0 )  £  V  if  and  only  if  there 
exists  a  nonnegative,  self  adjoint  operator  n  £  L(V) 
that  satisfied  the  algebraic  Rlccati  equation 


n  -  T*(J1  -  nB[R  *  B*I1B]  ^‘iDT  ♦  Q.  (3.2) 

If  for  any  control  sequence  that  makes  J(u)  finite 
the  state  x(k)  approaches  0  asymptotically,  then 
there  exists  at  most  one  nonnegative  self  adjoint 
solution  to  (3.2).  When  such  a  solution  to  (3.2) 
exists,  the  optimal  control  has  the  feedback  form 


u(k)  -  -F  x(k),  k  -  0.  1  .  2 .  (3.3) 


where  F  e  L(V,R  )  is  given  by 


F  •  [R  ♦  B*nB]_1B*nT, 


the  optimal  trajectory  satisfies 


u(k)  -  S  x(0),  k  -  0,  1,  2 .  (3.5) 


with  S  -  (T-BF)  e  L ( V )  and  -  <11  x(0).  x(0)>y. 

A  sufficient  condition  for  the  closed-loop  system 
to  be  uniformly  exponentially  stable  (i.e.,  for  the 
spectral  radius  of  S  to  be  less  than  one)  is  that  Q 
be  coercive  (i.e.,  Q  >  a  >  0). 


Since  F  e  L(V,  R  ).  there  exists 


t  -  (f.  f,  ...  f  )  e  Xv 

12  ”  J-l 


such  that  the  optimal  control  law  in  (3.3)  can  be 
written 


u(k)  =■  -<f.  x ( k )  >v,  k  -  0.  1,  2 .  (3.6) 

or  [u(k)]j  -  -<f.,  x(k)>V’  J  *  1.  2.  ...  m.  The 
vector  f  is  called  the  optimal  functional  feedback 
control  gain. 


The  optimal  LQG  state  estimator,  or  observer. 


x(k+l )  -  T  x(k)  +  B  u(k)  +  P  [y (k)  -  C  x(k)], 

k  -  0.  1.  2.  ...  <37) 

D 

Here,  x(0)  €  V  and  the  estimator  gain  F  e  L(R  ,  V) 


F  -  TnC* [R  ♦  cnc*]  1  (38) 

with  IT  in  L(V)  the  minimal  nonnegative  self  adjoint 
solution  (if  one  exists)  to  the  algebraic  Rlccatl 
equation 

n  =  T(fi  -  fic*[R  *  cnc*)"1cn)T*  *  q  , 

where  Q  e  L(V)  is  nonnegative  and  self  adjoint  and 
and  R  is  a  p  x  p  positive  definite  symmetric  matrix. 
If  Q  is  trace-class,  then  (3.7)  is  an  optimal 
Kalman-Bucy  filter  for  (2 . 10) -(2 . 11 )  if  a  zero-mean, 
stationary  white  noise  process  with  covariance 
operator  Q  is  added  to  the  right  side  of  (2.10)  and 
a  similar  noise  process  with  covariance  matrix  R  is 
added  to  the  right  side  of  (2.11)  (see  (1).  In 
compensator  design.  Q  and  R  often  are  chosen  to 
produce  certain  deterministic  properties  in  the 
closed-loop  system,  such  as  stability  margins  and 
robustness;  then,  Q  need  not  be  trace-class  Con¬ 
ditions  for  existence  and  uniqueness  of  solutions 
to  (3.9)  are  analogous  to  those  for  the  control 
Rlcattl  equation  (3.2). 


Since  F  e  L(R  ,  V),  we  can  write 


Fy  -  fT  y.  ye  RP. 


where  f  -  (f  l>  ...  f  )  €  X  V  is  known  as  the 

12  ■  j-i 

optimal  functional  estimator,  or  observer,  gain. 


The  optimal  compensator  consists  of  the  state 
estimator  in  (3.7)  and  the  control  law 


u(k)  -  -F  5t(k)  , 


k  -  0,  1.  2 .  (3.11) 


with  the  control  gain  F  given  by  (3.4)  with  the 
solution  to  (3.2)  and  the  estimator  gain  given  by 
(3.8)  with  the  solution  to  (3.9).  The  optimal 
closed-loop  system  then  consists  of  (2. 10)— (2.11) , 
(3.7)  and  (3.11).  It  can  be  shown  that  the 
spectrum  of  this  closed-loop  system  is  o(T  -  BF) 
o(T  -  FC) . 


Our  approximation  of  the  infinite  dimensional 
LQG  problem  and  compensator  is  based  on  approxima¬ 
tion  of  the  operators  T,  B,  C,  Q  and  Q  by  operators 
on  finite  dimensional  subspaces  VN  of  V.  We  denote 
by  the  orthogonal  projection  of  V  onto  V^,  and 
we  assume  that  PN  converges  strongly  to  the  iden¬ 
tity  of  operator  on  V.  Also,  we  assume  that  there 
exist  operators  TN,  QN>  QN  e  UVN)  and  QN  and  QN 
nonnegative  and  self  adjoint,  B^  e  L(R*,  V^)  and 
CN  c  L(Vn.  Rp)  such  that  yN.  T*Pj,.  QnPn  and  QfJPN 
converge  strongly  to  T,  T  ,  Q  and  Q,  respectively, 
and  Bjj  and  CjjPn  converge  in  norm  to  B  and  C, 
respectively  as  N  -  <». 


For  each  N,  we  define  a  finite  dimensional  LQG 
problem  for  the  open-loop  system  with  the  operators 
Tn.  BN  and  and  for  the  operators  Q^ ,  R  in  the 
optimal  control  problem  and  the  operators  Q^  and  R 
in  the  estimator  problem. 


The  solution  to  this  N  LQG  problem  requires 
the  solution  of  the  two  algebraic  Rlccatl  equations 
obtained  by  replacing  the  operators  in  (32)  and 
(3.9)  with  their  Nth  approximations: 


nN  =  T;,nN  -  nNBNlR*B‘nN8Nr  B‘VTN  *  QN  (312) 


"n  1  VnN  -  "nS^Wn1"  WTN  *  QN  ,313> 

The  optimal  control  and  estimator  gains  for  the 
problem  are  given  by 


fn  -  lR  -  b^)-1  B*nNTN 


h  ■  WV  WnJ’  <315) 

If  the  pairs  (T„,  B„)  and  (T*,  C*)  are  stabili- 

table.  (3.12)  and  (3.13)  have  nonnegative,  seif 

adjoint  solutions;  if  the  pairs  (Q^,  T^)  and 

(Q„ .  T*)  are  detectable,  the  nonnegative,  self 
N  N 

adjoint  solutions  to  (3.12)  and  (3.13)  are  unique 
and  positive  definite. 


w(t.O)  -  Dw(t.O)  -  0.  t  »  0.  (4.2) 

CQ  D2  w(t,l)  +E1  D2  w(t.l)  »  0,  t  i  0,  (4.3) 

cQ  D3  w(t,l)  ♦  El  D3  w ( t , 1 )  -  -u(t),  t  i  0,  (4.4) 
The  constants  in  (4.1)— (4.4)  are 

cQ  =  .0001,  El  -  .01333. 


As  in  the  infinite  dimensional  problem,  we  have 
-N 


functional  control  gains  f 


(f;  f2  ' ' '  f«)T  6  X  VN 

and  functional  estimator  gains  f  =  (f^  ?!*...  fN)T 

id  m 

e  X  corresponding  to  F^  and  F^,  respectively. 

We  will  give  the  formulas  for  computing  the  func¬ 
tional  control  and  estimator  gains  for  the  single 
Input/single  output  case  (m  =  p  -  1);  the  general 
Case  is  a  straightforward  extension.  When  m  *  1. 

B  and  B„  are  elements  of  V  and  V„,  respectively, 

N  N 

and  when  p  -  1,  C  and  CN  are  elements  of  V  and  V^, 
respectively.  The  following  two  formulas  then 
follow  from  (3.14)  and  (3.15): 


_N 


fl  ■  <R  *  ^W'1  Tn"nV 


(3.16) 


As  indicated  by  (4.4),  the  control  force  u(t)  is 
applied  to  the  right  end  of  the  beam.  For  the 
setup  in  Section  2,  the  appropriate  spaces  are: 

X  -  {<p  e  H2(0,1):  $(0)  >  D4><0)  -  0)  X  1^(0. 1)  (4.5) 

V  =  <$  e  H3  (0.1)  :  4>(  0)  =  D<J>(0 )  »  Dz4>(l)  *  0) 

(4.6) 

x(H1(0,l):  4>(0) )  -  0. 


W  -  (iJmJi)  e  X:  Cq<P  +  El <J>  e  H4(0,1). 

4>(0)  -  D<t> ( 0 )  =>  0,  cqD24»(  1  >  +  EXDZ«t>(  1 )  -  0) 
The  state  vector  x(t)  is 


(4.7) 


?N  *  ?i  -  <R  +  w S*'1  WnV  <3J7> 

For  the  example  in  Section  4,  it  can  be  shown 
that  the  closed-loop  operators  for  the  approxi¬ 
mating  LQG  problems  are  exponentially  stable, 
uniformly  in  N:  i.e.,  there  exists  constants  M  >  0 
and  r  <  1.  Independent  of  N.  such  that 

|(Tn  -  BNFN)k|  ♦  |(T*  -  FNCN)k|  <  Mrk,  (3.18) 

Vk  )  0 


x(t) 


w(t) 

w(t) 


and  the  operators  A.  T  and  r*  are 


(4.8) 


|-EI  D4  -  CQ  D4  J 

r :  W  -  R1,  r(4>.i)f)  -  -C0D34>(1)  -  EID3<M1). 

r*:  r1  -  w,  r+u  -  $0u. 


(4.9) 

(4.10) 

(4.11) 


This  condition  and  our  assumptions  about  the  strong 
convergence  of  the  basic  operators  T„,  T*  Q„,  Q.,, 
BN  and  CN  guarantee  (see  (3),  [4])  that  the  func¬ 
tional  control  gains  fN  and  f^  converge  in  V  to  f 
and  f.  respectively.  See  [ 3 ] ,  [4],  [5]  for  more 
details  on  convergence  analysis  and  numerical 
aspects  of  the  compensator  approximation  discussed 
here,  including  results  on  convergence  of  the 
approximating  compensators  and  closed-loop  systems. 

4.  EXAMPLE 

Consider  the  equation  of  motion  for  a  claaped- 
free  Euler-Bernoulll  beam  with  Kelvln-Volght 
viscoelastic  damping 

w(t.a)  ♦  cQ  D4  w(t.s)  ♦  El  D4  w( t , s )  -  0,  (4.1) 


where  is  the  third-degree  polynomial  that  satis¬ 
fies 

$0(0>  -  D$(  0 )  =>  D24>(1)  >  0,  D34>(1)  *  -1  /El  .  (4.12) 

Hence 

rr*:  r1  -  r1,  rr*  -  1.  (4.13) 

Now  let 

Dom(A)  -  N ( T ) ,  A  -  i|N(r)  (4.14) 

and  let  (T(t):  t  >  0)  be  the  semigroup  generated 
on  X  by  A.  Then,  as  in  (2.9), 

T  *  T(t)  B  •  (I  -  T(t))<J>0.  (4.15) 


0  <  s  <  1 .  t  >  0 


Note  that  R ( f )  c  N( A)  . 


.  For  nuaerical  solution  of  the  LQG  problea,  we 


used  two  scheaes  to  approxiaate  the  beaa.  The 
first  aethod  is  a  Galerkln  scheae  that  uses  cubic 
B-splines  as  basis  vectors  for  the  approxiaation 
subspaces  in  the  general  approxiaation  framework 
discussed  in  Section  3.  The  second  approach  is  a 
aodal  approxiaation,  which  uses  the  natural  aode 
shapes  of  undaaped  free  vibrations  of  the  beaa  as 
basis  vectors  in  approxiaating  the  infinite  diaen- 
slonal  control  problea.  We  have  not  carried  out 
all  of  the  details  of  the  convergence  proofs,  but 
it  appears  that  aethods  siailar  to  those  in  [5] 
will  yield  the  required  proofs. 

For  the  optiaal  linear-quadratic  regulator 
problea  we  take  T  «  .01.  R  *  .01  and  Q  such  that 

<Qx(t) .  x ( t ) >v  -  |x(t)|x2.  (4.15) 

The  functional  control  gain  f  in  (3.6)  has  the  fora 
f  *  (fa.  f2)  e  V.  (4.16) 

In  the  spline-based  scheae,  the  function  4>Q  is 
an  eleaent  of  each  V^,  so  that  essentially  we  need 
only  approxiaate  the  sealgroup  T(t)  using  the 
Galerkln  aethod.  For  each  N.  we  solve  the  finite 
diaenslonal  Rlccatl  equation  (3.12)  by  standard 
nuaerical  aethods  and  coapute  the  Nth  approxiaation 
to  f  according  to  (3.16).  Figures  1-3  show  our 
approxlaations  to  f ^ ,  f”  and  f2  for  N  =  2,  4,  6,  8 
eleaents . 

Soae  of  the  splines  used  in  this  scheme  are  in 
X  but  not  in  V,  so  we  are  not  guaranteed  V- 
convergence  for  the  approxiaating  functional  gains. 
Figures  1-3  demonstrate  that  the  approxlaations  to 
f  converge  in  X.  It  seems  unlikely  that  these 

approxlaations  converge  in  V;  V-convergence  would 
N" 

require  that  f^  (1)  converge  to  zero,  but  it  does 
not. 

In  the  second  approxiaation  scheae,  the  natural 

aode  shapes,  which  are  used  as  basis  vectors  for 

the  spaces  V„  in  Section  3,  are  eleaents  of  V,  so 
”  N 

we  should  get  the  functions  f  converging  in  V. 

It  is  straightforward  to  coapute  the  operators  T^, 


Bjj  and  Qn  for  the  finite  diaenslonal  Riccati 
equations.  Figures  4  and  5  show  the  approxiaating 
functional  control  gains  for  N  *  10,  14,  18,  20 
aodes,  and  the  pointwise  convergence  of  fjN  is 
consistent  with  V-convergence  of  the  approxlaations 
to  the  functional  control  gain  f. 


Figure  4.  Modal  Approximations  to  fj 
N  *  10,  14,  18.  20  Modes 


Figure  6.  Modal  Approximations  to  li 

N  »  14,  18.  20,  22.  24  Modes 


Figure  5.  Modal  Approximations  to  I2 
N  -  10.  14.  18.  20  Modes 


Figure  7.  Modal  Approximations  to  fj 

N  »  14,  18.  20.  22,  24  Modes 


For  state  estimation,  we  take  the  measurement 


y ( t )  >  w(t.  .5); 


l.e.,  we  place  a  strain  guage  at  the  midpoint  of 
the  beam.  We  need  the  state  vector  to  be  in  the 
space  V  In  (4.6)  so  that  a  polntwise  measurement  of 
the  -second  derivative  of  the  displacement  will  be  a 
continuous  measurement  of  the  state.  In  the  esti¬ 
mation  problem,  we  take  R  -  1  and  Q  -  1 


We  use  only  the  modal  approximation  to  solve 
for  the  estimator  gains.  For  each  N  (number  of 
modes),  we  solve  the  finite  dimensional  estimator 
Rlccatl  equation  (3.9)  and  then  compute  the  approx¬ 
imating  functional  estimator  gains  f **  »  (fN,f|!) 

1  N-  z 

according  to  (3.17).  Figures  4-6  show  f"  and 

N  * 

fj,.  As  expected,  these  functions  converge  uni¬ 
formly  polntwise.  This  Is  consistent  with  V- 

n 

convergence  of  f  . 
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ABSTRACT 


'This  paper  presents  a  new  derivation  of  a  least-squares  lattice  filter  and  applies  the 
filter  to  identification  of  flexible  structures.  'Die  vector-channel  lattice  here  --  derived  in 
an  infinite  dimensional  history  space,  without  matrix  manipulations  or  geometric 
arguments  —  can  constrain  the  AR  coefficients  for  several  outputs  to  be  the  same. 
Numerical  results  for  a  simulated  flexible  structure  compare  the  vector-channel  lattice 
with  the  standard  lattice  filter.  These  results  show  that  the  frequencies  and  damping 
ratios  for  the  most  significant  modes  can  be  identified  adaptively  with  lattice  filters. 
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1.  INTRODUCTION 


CH, 

.35 

;SS 


LJB 

cvq 

#] 

rJ 


iJn 


•‘a* 

•v 

■:5 


In  recent  years,  adaptive  lattice  filters  have  been  used  increasingly  in  signal  processing  to  implement 
least-squares  algorithms  and  other  estimation  methods.  One  advantage  of  a  least-squares  lattice  filter  is 
that  it  is  recursive  in  both  time  and  order,  while  the  classical  least-squares  algorithm  is  recursive  in  time 
only.  Also,  lattice  filters  have  been  shown  to  possess  desirable  numerical  stability  [1],  Lee.  Morf  and 
Friedlander  [2,3]  developed  the  time  update  equations  needed  for  adaptive  lattice  filters  and  studied 
applications  to  speech  processing  and  adaptive  control.  For  derivations  and  discussions  of  lattice  filters, 
see  [4, 5, 6, 7], 

Potentially  important  new  applications  of  lattice  filters  appear  to  lie  in  adaptive  identification  and 
control  of  large  flexible  structures.  Because  such  sturcturcs  have  many,  usually  infinitely  many  modes  of 
vibration,  the  problem  of  determining  the  order  of  a  model  to  use  for  on-line  identification  and  control  is 
a  major  obstacle.  A  lattice  filter,  with  its  order-recursive  property,  can  in  principle  identify  the  effective 
order  of  the  structure  and  adaptively  increase  this  order  to  accomodate  new  excitations  to  the  structure  or 
decrease  the  order  as  faster  transients  in  the  structure  arc  damped  out  Sundararajan  and  Montgomery 
[8,9]  and  NViberg  [10,11]  were  fust  to  apply  lattice  filters  to  identification  of  flexible  structures. 
Sundararajan  and  Montgomery  used  the  lattice  to  estimate  the  number  of  excited  modes  of  a  flexible 
structure,  while  using  fast  Fourier  transforms  to  identify  the  frequencies  and  damping  ratios.  Wiberg  has 
studied  identification  of  structural  frequencies  and  damping  with  lattice  filters.  For  structural  identification 
from  free-response  data,  Wiberg  has  developed  a  "vibration  lattice",  which  enforces  the  constraint  that 
each  scalar  measurement  channel  for  the  same  structure  have  the  same  auto  regressive  (AR)  coefficients. 

To  understand  the  motivation  for  Wiberg's  vibration  lattice  and  our  vector-channel  lattice,  consider 
the  usual  Unear  AR  (auto  regressive)  model 

n 

y(>)=  ^  Ajy(l-j)  +  £„(0,  (1.1) 

/=! 

where  y(t)  is  an  m- vector,  each  Aj  is  an  m  x  m  matrix  of  parameters  to  be  estimated  and  £„(/)  is  the  error. 
Die  least-squares  lattice  filter  yields  a  least-squares  estimate  for  the  Aj 's  that  is  recursive  in  both  t  and  n. 
Now  suppose  that  y(i)  represents  m  sensors  on  the  same  flexible  structure  and  that  N  observable  modes 
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of  the  structure  are  excited.  Then  the  AR  mode!  above  with  n  something  less  than  2N  should  fit  the  data. 
However,  only  in  special  circumstances,  will  the  order  of  the  minimum  AR  model  that  fits  the  data  equal 
the  order  of  the  structural  vibrations,  since  the  AR  model  has  nm  eigenvalues  while  the  structural 
response  has  2N  eigenvalues.  In  general,  the  minimal  AR  model  will  have  extraneous  eigenvalues.  If 
modal  frequencies  and  damping  ratios  are  desired,  their  identification  will  be  complicated  severely  by  the 
extraneous  eigenvalues  of  the  AR  model. 

On  the  other  hand,  if  the  Aj 's  in  the  ARMA  model  are  constrained  to  be  scalars,  then  the  minimum 
n  is  always  2N.  This  constraint  represents  significant  information  about  the  system;  it  amounts  to  wnting 
the  scalar  AR  model  for  each  component  of  the  measurement  vector  and  imposing  the  constraint  that  all 
m  of  these  processes  have  the  same  AR  coefficients  because  they  are  measurements  of  the  same  linear 
system.  When  the  measurements  are  noisy,  incorporating  this  information  into  the  estimation  algorithm 
is  even  more  important  than  resolving  the  relation  between  the  number  of  modes  and  the  order  of  the  AR 
model.  Using  all  information  about  the  system  is  especially  important  in  on-line  identification  of 
distributed  systems  like  flexible  structures.  Even  if  the  sensor  data  is  virtually  free  of  random  noise,  the 
number  of  excited  modes  often  is  greater  than  the  allowable  order  of  the  lattice  filter,  so  that  the  higher 
frequency  modes  amount  to  unmodeled  distrubance  in  the  sensor  data. 

The  lattice  derived  here  solves  the  problem  of  constraining  the  AR  coefficients  for  multiple  outputs 
to  be  the  same.  We  refer  to  this  lattice  as  a  vector-channel  lattice  because  it  treats  several  outputs  as 
components  of  one  vector  measurement  channel,  as  opposed  to  previous  lattices,  which  treat  each  scalar 
output  as  a  single  channel.  In  general,  there  can  be  multiple  vector  measurement  channels.  When  each 
channel  contains  only  one  component,  the  vector-channel  lattice  reduces  to  the  standard 
multi-input-multi-output  lattice,  which  we  refer  to  as  either  the  standard  lattice  or  the  scalar-channel 
lattice. 

Although  application  of  the  classical  recursive  least-squares  algorithm  to  the  case  where  the 
measurements  are  vectors  but  the  AR  coefficients  are  scalars  is  straight  forward  and  well  known  (see  [7]), 
a  lattice  filter  for  this  problem  apparently  was  not  derived  until  Wiberg  introduced  his  vibration  lattice  for 
structural  identification.  The  vector-channel  lattice  differs  from  the  vibration  lattice  in  two  respects,  f  irst, 
one  coefficient  that  is  a  scalar  in  Wiberg 's  and  previous  lattices  is  a  matrix  in  the  vector-channel  lattice, 
the  dimension  of  this  matrix  being  determined  by  the  dimension  of  he  vector  measurement  channels. 
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Second,  the  vibration  lattice  docs  not  accomodate  an  input,  or  control;  i.c.,  the  moving  average  part  of  an 
AR\1A  model.  This  capability  is,  of  course,  important  in  adaptive  identification!  and  control.  As  in 
previous  lattices,  a  moving  average  must  be  treated  as  an  additional  measurement  channel  in  the 
vector-channel  lattice.  The  trick  is  to  constrain  the  AR  coefficients  to  be  the  same  for  all  true 
measurement  channels  while  allowing  different  coefficients  for  the  input  channels.  Our  method  for 
including  the  MA  terms  in  the  vector-channel  lattice  is  discussed  in  Section  3.3. 


In  Section  2,  we  derive  the  vector-channel  lattice  filter  and  the  algorithm  for  generating  the  AR 
coefficients  from  the  lattice  data.  This  derivation  differs  from  previous  derivations  of  least-squares  lattices 
in  several  importants  respects,  all  of  which  stem  from  the  fact  that  we  work  in  a  Hilbert  space  of  infinitely 
long  history'  vectors.  Previous  lattice  derivations  have  dealt  with  finite  history  vectors  of  changing  length 
and  have  relied  heavily  on  matrix  manipulation  and/or  geometric  arguments,  including  nonorthogonal 
projections.  We  use  neither  type  of  argument;  our  derivation  requires  only  elementary  properties  of 
orthogonal  projections  in  Hilbert  space. 


The  original  motivation  for  the  derivation  here  was  not  to  develop  a  vector-channel  lattice.  In 
[12,13]  we  developed  infinite  dimensional  AR\1A  models  for  certain  classes  of  distributed  systems,  with 
one  aim  being  to  use  these  ARMA  models  as  the  basis  for  an  approximation  theory  that  will  predict  the 
behavior  of  lattices  of  increasing  finite  order  in  adaptive  identification  and  control  of  infinite  dimensional 
systems.  In  such  analysis  (which  is  continuing),  infinite  histories  of  inputs  and  outputs  are  used  Hence, 
the  infinite  dimensional  history  space  in  which  we  work. 


In  Section  3,  we  summarize  the  vector-channel  lattice  algorithm  and  discuss  its  implementation. 
From  the  derivation  of  the  lattice,  we  deduce  an  algorithm  for  one-step-ahead  prediction  and  show  how¬ 
to  make  copies  of  an  input  so  that  in  a  vector  ARMA  process  the  AR  coefficients  can  be  constrained  to 
be  the  same  for  each  channel  but  not  the  the  MA  coefficients. 


In  Section  4,  we  compare  numerical  results  obtained  with  standard  and  vector-channel  lattices  for  the 
simulated  forced  vibrations  of  a  flexible  structure  with  two  position  measurements.  First,  we  use  the 
standard  lattice  to  treat  the  measuremets  as  two  independent  ARMA  processes.  Then  we  use  the 
vector-channel  lattice  to  treat  the  measurements  as  the  two  components  of  one  vector  process,  thereby- 
imposing  the  constraint  that  each  measurement  process  have  the  same  AR  coefficients.  We  use  the  lattice 
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to  identify  both  the  number  of  significantly  excited  modes  and  the  frequencies  and  damping  ratios  of  those 
modes.  We  give  results  for  simulations  both  with  and  without  noisy  measurements. 
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2.  DFRIVA  TION 


Definitions  and  Order  Updates 


The  history  space  in  this  paper  is  the  following  Hilbert  space  of  infinite  column  vectors: 


T  T  T  T 

=  [  xt  x2  x2  ...  ]  reach  X:  is  a  real  m-  vector 


and  ||i/' ||  =  <  \p  ,ip  >  <  oo 


where 


OO 

<  'I'.'i'  >  =  ]T  ^_l  xfxj 


(2.1) 


J=\ 


and  A  (the  forgetting  factor)  is  a  positive  real  number.  Throughout  this  paper,  <  .,.  >  means  the  inner 
product  in  (2.1). 


1  he  following  derivation  assumes  a  sampled  process  with  p  channels,  each  of  which  contains  m  scalar 
measurements.  In  other  words,  the  measurement  from  channel  i  at  time  t  is  the  m-vector 


/(<)=[  .I’ju)  ...V>>  f.  1=  1  (2.2) 

The  inlimte  history  vector  of  channel  i  is 

=  [  (y‘U))T  ( y‘(t  -  1)  )r  .  f  .  (2.3) 


In  applications,  all  hut  a  finite  number  of  terms  in  ip\t)  can  be  set  to  zero.  The  main  advantage  of  the 
infinite  dimensional  history  space  for  the  lattice  derivation  is  not  that  it  allows  histories  of  infinite  length, 
but  that  it  accomodates  histories  of  all  finite  lengths  in  the  same  space. 

ITie  following  definitions  are  necessary.  For  any  integer  t  and  nonnegative  integer  n: 

S(<)  =  span  {  i/i'u),  ....  i//p(t)  }  (2.4) 

//(, (/)={<>}.  //„(!)=  S(1 2)  ®..®S{t-n)  (2.5) 

Pn(t)  =  Orthogonal  Projection  onto  //„(/)  (2  6) 


-  5  - 


'A-y‘jw(.vV 


av, -v 


v.v 


Eq(i)=S(i),  E?nd)  =  [  /  -  Pn(t)]  S (/).  forward  error  space 


(2.7) 


E$(t)=S(i),  £*(/)=[/  -  /*„(/+  1)]  S{i  -  n) ,  fackward  error  space  (2.8) 

P{(t)=  Orthogonal  projection  onto  E}n(i) ,  Pb( t)  =  Orthogonal  projection  onto  Eb({)  (2.9) 

Note  that  the  definition  of  the  backward  error  space  E%(t  -  1)  amounts  to  a  Gram-Schmidt 
orthogonalization  of  the  history  spaces  S(t-l),  ...  ,  S(t-n). 

The  rest  of  this  section  derives  a  number  of  recursions  from  the  definitions  given  so  far.  Like 
(2.4)-(2.9),  subsequent  equations  in  this  section  are  valid  for  any  integer  t  and  any  nonnegative  integer  n, 
except  where  otherwise  indicated.  Of  course,  for  implementation  of  the  lattice  algorithm,  initialization  at 
some  finite  initial  time  is  necessary,  and  this  is  discussed  in  Section  3.  In  the  current  section,  initialization 
need  not  be  considered. 

A  first  consequence  of  the  foregoing  definitions  is 

//„+,(')=  HnU~  <)©  l&t-  1  )=//„</)©  /.«V  -  1)  (2.10) 

l 

where  the  notation  ©  indicates  the  direct  sum  of  orthogonal  subspaccs  of  £2(/vm,A).  Because  of  (2.10) 


and  since  Ef(l)  and  T*(t  —  1)  arc  orthogonal  to  ff„d)  ,  elementary  properties  of  projection  yield 

p„+i  U)=P„U-  1)  +  liu-  1)  =  /y  n  +  pb„d-  1)  (2.11) 

Pid)  Pn(t)=Pbn(t-  1)  />„(/)=<).  (2.12) 

From  (2.1 1 )  and  (2. 12),  it  follows  that 

[  /  -  Pn+](t+  !)]=[/-  fid)  ]  [  /  -  Pnd)  ],  (2.13) 

[  /  -  /WO  ]  =  [  /  -  Pbnd- I)  ]  [  /  -  Pnd)  ]  (2.14) 

The  definitions  of  the  forward  and  backward  error  vectors  are  similar  to  (2.7)  and  (2.8): 

f‘d.t)='l'X0  .  /‘n(0=  [  /  ~  Pn(t)]  'I'X'),  Forward  error  (2.15) 

b^t)  =  \p‘(t)  ,  b‘nd)=[l  -  P„(t  +  l)]t l>\l  -  n) ,  i=  1,  ...p  Backward  error  (2.16) 


w 


E]n(i)  =  span  {fnU) }  and  E„(t)  =  span  {  b‘n{t) }  ,  i  =  1,  ...p. 


The  order  update  formulas  for  the  forward  and  backward  error  vectors  follow  from  (2.13)-(2  16) 
Substituting  (2.14)  into  (2.15)  for  n+  1  yields 

/U»(0=[/  -  />„+,(')]  ^(0=[/  -  Pn('~  D][/-  Pn(‘)  ]  <A'(0, 


which  in  view  of  (2.7)  becomes 


/Ui(o=  [/  -  -  /*&-  D/n(o 


for  i  =  1,2,  ...  p.  Similarly,  (2.13)  and  (2.16)  with  (2.8)  yield 

ClW=K('-  l>  -  KiObnU-  1). 


To  make  these  order  update  equations  useful,  the  following  elementary  and  standard  result  is  required 
(see  [14],  page  56). 

lemma  2.1.  Let  II  be  a  real  Hilbert  space,  E—span{h],  h2 . hp  }  a  subspace  of  II  and  P£  the 

orthogonal  projection  of  II  onto  E  .  Then  for  any  h  e  II, 


Pi;  ft  =[  ft i  h2  ■■■  ftp  ]  a  (  =  V  al  h:  ) 


where  a  e  Rp  satisfies 


Rf  a=  d 


with  /?£  the  nonnegative  symmetric  matrix  whose  (i,j)  clement  is  <  hb  hj  > u  and  d  the  p-vector  whose  ith 
element  is  <  hb  h>  [{  .  For  each  he  II,  at  least  one  solution  to  (2.21)  exists.  The  matrix  RE  is  positive 
definite  if  and  only  if  the  vectors  A,  are  linearly  independent. 

The  following  matrices  are  useful  for  applying  Lemma  2.1  in  (2. 18)  and  (2. 19).  For  n  =  0,  1,  ... 


Kn+\(t)  =  p  x  p  matrix  whose  (i.j)  element  is  <  f‘n{() ,  bJn(t~  I)  >  , 
Ren{i)  —  p*  p  matrix  whose  (i.j)  element  is  <  >  , 


Rrn(t)  =  p  x  p  matrix  whose  (i.j)  clement  is  <  b'n(t) ,  b„(t)  >  . 


(2.24) 


Also,  we  define 

[  /„(')  ]  =  [  /i( ')  flu)  fPn  ].  [  bnU)  3  =  [  b\{l)  b2n{t)  ...  %  ]  ;  (2.25) 

i.e.,  ]  is  the  matrix  containing  the  p  infinitely  long  columns  /'„(()  ,  and  [  bn(l)  ]  is  the  corresponding 

matrix  containing  the  bn(t)'s.  With  these  definitions  and  (2.17),  Lemma  2.1  implies  that  (2.18)  and  (2.19) 
are  equivalent  to 

[  /*+,(')  3  =  C  /«(')  3  -  C  bn(t~  1)  )  R-\t~  1)0')-  (2.26) 

[  bn+l(t)  ]  =  [  bn(t-  1)  ]  -  [  fn(t)  ]  R-n\t)Kn+x(t).  (2.27) 

Here,  KnJrX(l)  means  any  pxp  matrix  a  such  that  R*(t)  a  =  Kn+  j(r)  and  R~^(t-  l)A'„r+l(t)  means 

any  pxp  matrix  /?  such  that  Rrn(t  —  1)  (S  =  Kj[+](t).  That  at  least  one  such  a  and  one  such  /(  exist  follows 

from  Lemma  2.1  and  the  existence  of  the  projections  in  (2.18)  and  (2.19). 

While  (2.26)  and  (2.27)  are  more  explicit  than  (2.18)  and  (2.19),  they  still  are  not  useful  for 
computation  because  columns  of  the  matrices  [/j,(f)j  and  [  bn(i)  ]  are  infinitely  long.  It  turns  out  that 
only  the  top  m  rows  of  (2.26)  and  (2.27)  are  used  in  the  lattice  filter.  Hence,  it  is  useful  to  define  the 
following  mxp  matrices: 

e„(f)=  Top  m  rows  of  [  fn(l)  ],  r„(/)  =  Top  m  rows  of  [  bn(t)  ]  (2.28) 

Then,  the  top  m  rows  of  (2.26)  and  (2.27)  arc 

e„+,(0=  en(t)  -  rn(t-  \)  R~r(t  -  forward  residual  (2.29) 

rn+1(t)=  r„(t-  1)  -  en(r)/?“tf(r)A"n+|(t),  backward  residual.  (2.30) 

These  two  equations  are  used  in  the  lattice  filter  for  order  updates  of  the  forward  and  backward  residuals. 
For  order  updates  of  R^i)  and  Rr(t),  substituting  (2.29)  and  (2.30)  into  (2.23)  and  (2.24)  yields 


*;.,(/)= -  k„+.(')/C<'-  nOo, 


C)(0=  RrnU-  1)  -  Kl+](t) R'e{t) Rn+M 


l  ime  Updates 


The  following  vectors  are  elements  of  f2(^m>  '  )  '■ 


position  1 


0'  =  [  0  0  ...  1  0  ...  0  ..  ]r,  /=  1.2 . 

4>‘nU+  D=[/-  P„0  +  1)  ]  <t>‘  .  /  =  1.2,  ....  m; 


where  Pn(l  4-  1)  is  the  projection  onto  Hn(t  4-  1)  ,  as  before.  The  subspacc 


Hn{t  4  1)  =  span  { \p  ’(/),  ....  i pp(t),  ....  i p\t  -  n  +  1) .  ipp(t  -  n  +  1) . <?m  J 


>•  < 


is  useful  in  deriving  the  time  update  equations  for  the  variables  Kn(l)  .  I  roin  (2.3)-(2.5)  it  follows  that 


//„(  t  4-  1 )  =  span 


(  0(m)  )  (  0(m)  ) 
\  +  \t-  1)/’  -  \+pU-  1) )' 


Ip  (I  -  n)  r  1  i pp{l  -  n) 


Hn(l  +  0=  +  1)©  span  {  <j>n(t  +  1),  ....  <p”(i  4-  1) }  , 


where  is  the  mxm  zero  matrix.  If 


Pn(t  4-  1)=  Projection  onto  Hn(t  4-  1)  and 


Qn(t  +  1)=  Projection  onto  span  {  <pn(t  4-  1),  ...,  4>™{t  4-  1)} 


it  follows  from  (2.34)  that  Pn(t  4-  l)  =  Pn(t  4-  1)  4-  Qn(t  4-  i )  ,  so  that 


[/-  Pn(t  +  !)]  =  [/  -  Pn(l  +  1)]  +  Q„(t+  1). 


Tr^r.’Tf.  v»  v7V-.--.r.-v\-VYsr_  VV  w  *V vvv-.-  VV v.  IT  .-V. 


Substituting  (2.36)  into  (2.16)  yields 


/>'„(')=  n  +  r 


where 


II  =  [  /  -  />*(/+  D]  *J(t-n),  L  =  0„(r  +  1)  *'(*-«). 


From  (2.3)  and  (2.31)  it  follows  that 


1 1/(1-  n)  = 


t//(f-  n-  1) 


rn 

+  ^  n)4>k  . 


As  in  (2.2),  is  the  k,h  sc;ilar  measurement  in  channel  j.  Also,  in  view  of  (2.33), 


H  =  [  /  -  Pn(t  +  D] 


—  n—  1) 


=  /  0im>  \  /  0'm' 
VC/-  WltH'-*-  I)/  \^(f-  ') 


According  to  (2.35)  and  Lemma  2.1, 


where 


r  =  [  <£„('  +  i)  ]  g-W, 


[  </>„«  +  l)  ]  =  [  0i(f+  1)  4>l(t+  1)  ...  K(‘  +  l)  ] 


y  A 

G„(t)  =  mx  m  matrix  whose  (i,j)  element  is  <  <£'„(/+  1),  <pJ„U+  1)  > 


d  =  m- vector  whose  k  clement  is  <  4>„(t+  1),  >  . 


The  expression  G~ ’(/)  d  means  any  m-vector  a  such  that  Gn( t) a  =  d  ,  and  it  follows  from  Lemma  2.1  that 
at  least  one  such  a  exists.  Since  [  /  —  Pn(t  +  1)  ]  is  self-adjoint  (being  an  orthogonal  projection), 


k  element  of  d  =  <  4>k(i  -I-  1),  \p  \l  —  n)  > 

=  <  <f>kl  l  -  P„(t  +  1)  ]  *  \l  -/.)>  =  <  4>k  ,  bJn(t)  >  =  rkJ(t)  ■ 


i.e.,  the  vector  d  in  (2.42)  is  the  jth  column  of  the  matrix  rn(t). 


Now,  since  [  /  —  Pn(t  +  1)]  is  an  orthogonal  projection,  it  is  self-adjoint  and  equal  to  its  square. 


Hence,  from  (2.37), 


</„(!+  D,  bjn(i)>  =  <ip\t+  i),  bjn(t)>  =  <^u+  i),  n>  +  <^‘u+  i) .  r >  .  (2.44) 


According  to  (2.1)  and  (2.38), 

<  ^(i+  i),  n  >  =  2  bid-  i)  >  =  ;.  <f‘n(t ),  bid-  i)  > 


<  ip\t+  1),  4>knd+  i)  >  =<  [/-  /yr  +  0]</A'  +  1) ,  <t>k  > 

=  </'„«  +  i).  4>k  >  =  4V+  i). 

using  I  from  (2.39)-(2.43)  yields 


<  \jj\t  4-  1) ,  P  >  =  [  i,h  column  of  e„(f  +  1)  ]r  Gn  V)  [/*  column  of  rn(t)  ] 


Now,  equations  (2.44)-(2.46)  show  that  the  time  update  equation  for  A.'„+|(0  (defined  by  (2.22)  )  is 


K„+\U+  D=  AKn+]d)+  e^d+  l)Gn\t)r„d). 


The  final  equation  needed  to  complete  the  residual  error  lattice  is  the  order  update  equation  for 
<7„(r)  f  rom  (2.32)  and  (2.41),  the  elements  of  this  matrix  are 

G‘nJd)=<  4>i(i+  l),  &U+  l)  >  =  <  <4('  +  1).  >  (247) 

=  <  <p‘ ,  >  -  <  Pnd  +  l )  <t>‘ .  >  • 


The  subspace  lln(t  4-  1)  can  be  written  as  the  following  direct  sum  of  orthogonal  subspaces: 

nnd  +  l)=C,(oiil2(0©  -  ©4(0, 


so  that 


pn(t+  d=4-i(o  +  4-2(0  +  -  +  4(0 , 


According  to  (2.17),  (2,28)  and  Lemma  2.1, 


Since  ]  is  the  error  remaining  after  orthogonal  projection  of  the  data  taken  through  time  t  (recall 
(2.2)-(2.3))  onto  the  history  space  //„(t),  the  coefficient  matrices  AnJt)  rmninve  the  norm  of 

[/,(!)]  over  all  autoregressive  models  of  order  n.  ITie  matrices  Anj^t)  are  called  the  AR  coefficients  at 
time  t  . 

Substituting  (2.52)  and  (2.53)  into  the  order  updates  in  (2.26)  and  (2.27)  and  matching  coefficients  of 
the  history  vectors  4*\t  —  j)  yields 


'K+\j(t)=AnJ(t)-BnJ(t-  1  )R-r(t-  1  )<+I(0. 

j  =  1,2 . n,  n  =  1,2,... 

(2.55) 

An+Un^(t)=R~\t-  1  )*£,</). 

n=  0,1,2. ... 

(2.56) 

Bn+ 1  j+ ,  (/)  =  BnJ(t  -  1 )  -  Anj  0)  R-\t)  Kn+  ,(r) 

j  =  l,2,...,n,  n=  1,2,... 

(2.57) 

l*n+uU)=  R;eU)Kn+](t), 

=  0,1,2, ... . 

(2.58) 

The  AR  coefficients,  An+](l)  ,  can  be  generated  with  (2.55)-(2. 58);  however,  this  algorithm  requires 
that  the  AR  coefficients  be  calculated  at  every  sampling  time  because  Bnl(l  —  1)  is  needed  to  compute 
-4„+|(t).  The  following  derivation  provides  an  algorithm  for  computing  the  AR  coefficients  at  any  time  l 
without  the  values  of  Anj(t  —  1)  and  BnAl —  1)  . 

Using  (2.1 1)  in  (2.32)  with  the  notation  in  (2.40)  yields 

4>n+\U+  !)=[/“  Pn+i«+  1  )]<*>'  =  [/-  Pn('+  l)]</>''  -  Pbn(D4>‘ 

rn\‘) 

=  fat  +  1)  -  =  fat  +  1)  -  t  bn(t)  ]  Kr(‘) 


[<W'  +  1)  ]  =  C  <*>„('+  D]  - 


Since  Pn(t  +  \)<t>‘  e  //„(/  +  1)  ,  it  follows  from  (2.5)  and  (2.32)  that  there  exist  p  x  m  matrices  CnJi) 
such  that 


.  r  . 


[<*>*+>('  +  1)]  =  0  ]  -  y  [>/'('-y+  D]  C„j(‘),  n=  1,2,... 


where  [<£]=[</>'  <£2  ...  4>m  ]  .  Substituting  (2.53)  and  (2.60)  into  (2.59)  and  matching  coefficients  of 


the  history  vectors  yields 


C„^j(‘)=Cnj(i)  -  BnJU)Rn\t)rnr(t)  /=1.2 . «=I,2 . 


,.*+,(')  =  O'KV).  n=  0.1.2 . 


The  expression  [  bn(t)  ]  will  denote  the  part  of  [  bn(t)  ]  below  the  first  m  rows.  Then,  from  (2.53) 


[£,(*)]  =  W-«-  1)]  -  y  [^(t-j)-]BnJ(t),  n=  1.2 . 


Next,  the  developcmcnt  in  (2.37)-(2.43)  yields 


MO  ]  = 


[6ji_  u]  +  0*('+  i)]o;'(0M0. 


After  substitution  of  (2.60)  into  (2.64),  the  new  equation  from  row  (m+  1)  down  is 


[  MO]  =  [bn(i-  1)]  -  y  [Mt-/)lCn.U)O;\t)rJ0.  n  =  1,2, 


Finally,  equating  the  right  sides  of  (2.63)  and  (2.65),  using  the  right  side  of  (2.53 )  for  bn(t  -  I),  and 


matching  the  coefficients  of  the  history  vectors  yields 


l)=Bnj(t)  -  Cnj(t)  G„  \()rn{t) ,  j  =  1.2,....  n,  #i  =  1.2 . 


I  he  algorithm  consisting  of  equations  ( 2. 55)-( 2. 58),  (2.61),  (2.62)  and  (2.66)  can  be  started  at  anv 


time  to  obtain  the  AR  coefficients.  The  order  of  the  AR  coefficients  need  not  remain  the  same  from  one 


time  to  another.  In  addition  to  this  flexibility,  this  method  is  computationally  attractive  when  the 


sampling  rate  is  high  but  the  plant  and/or  control  parameters  need  not  be  updated  at  every  sample  time. 
The  residua]  filter  is  the  only  part  that  must  be  computed  at  every  time  step. 
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3  IMl’I  BlliM'.VI  ION 

Summary  ot  tlic  1  atncc  Alcorithms 

Recall  the  definition  (2.2)  for  the  measurement  vectors  .  We  dclinc  the  m  ■  p  matrix 

y(/)  =  [vV)  v-2(0  -  /«)] 

which  contains  the  measurements  at  time  t.  We  start  the  lattice  tiller  at  time  t  -  0  with  the  following 
initialization: 

R^(  —  1 )  =  A’,)(  —  1 )  =  <x  /  ,  A„  „,«)=<>  for  w  +  I  >  «  >  0  .  (3.1) 

where  u  is  a  sm;dl  number.  If  a  =  0,  this  initialization  corresponds  to  what  commonly  is  called 
prewindowing,  which  means  taking  y(t)  =  0  for  all  t  <  0.  We  have  found  that  taking  a  on  the  order  ol 
If)-7  avoids  numencal  problems  with  inverting  certain  matrices  in  the  filter,  as  discussed  below. 

Residual  hrror  1  attice 

for  each  t  >  0: 


e0(r)  =  r„ (/)  =  ;■( /) .  A,,',U)=  A’^t / )  =  dn/  i\t)  +-  z  R^t  -  1), 

<'u(0  =  ’ 

(3.2) 

for  n  =  0  to  t- 1 : 

jfntl(/)  =  ;.A',t ,(/-  i)  f  <•>)<;;'( t-  i ku-  n 

(3.3) 

e„+1(f)=  en(t)  -  rn(i  —  \)RnrU~  I)  A„7+,(r) 

(3.4) 

rn+iU)  =  rnU-  I)  -  <’,(/)  /?;'(/)  An+,(r) 

(3.5) 

<+1(')=<(0  -  A,+  l(f)R;,(/-  DA„;+1(/) 

(3.6) 

<+,(/>=<('-  1)  -  ^^(OCW^oW 

(3.7) 

«„+.(')=  <6,0)  -  rn(0An“V)rJ(t)  #.=  0,1,2... 

(3.8) 

I  he  maximum  n  can  be  increased  by  one  with  each  successive  time  step,  in  practice,  though,  n  is 
not  increased  past  some  fixed  order,  due  to  finite  computational  capacity. 


AR  Coefficients 


The  top  m  rows  of  (2.52)  are 

n 

c„(t)  =  y(t)  -  V  y(t~j)Anj(t)  .  n=  1,2,  ... 
i=  i 

At  any  time  t,  the  AR  coefficients  Anj(t)  can  be  generated  from  the  residual  filter  data  by  the  following: 
lorn  =  1 ,  t- 1  andfori=!.n 

Cn+U(t)=  CnJ{t)  -  Hn  i(t)  R~\t)  rj(l)  (3.9) 

«„+!;+,  (0=  [  «*.,(')  ~  ^(OC-VO^O]  -  Ani(t)R;eU)Kn+i(t)  (3.10) 

An+uU)=A „,(<)  -  [«„,(/)  -  Cni(t)C,A(i)rn(i)l  R?U  -  D  *„+,(»)  (3.1  D 

where  equations  (2.56),  ( 2. 5S),  and  (2.62)  arc  used  as  end  conditions. 

Matrix  Inverses 

It  A’£(t),  Rr„(i)  and  (/n(i)  are  nonsingular,  then  R„e(i),  R^r(t)  and  G^\t)  are  the  respective  inverses.  If 
any  of  these  matrices  is  singular,  then  the  inverse  should  be  interpreted  as  indicated  alter  (2.27).  If  the 
measurements  are  independent,  then  R%t).  Rrn(i)  and  Gn(i)  are  positive  definite  unless  one  of  the  channels 
(or  a  linear  combination  of  the  channels)  fits  the  nlh  order  AR  model  exactly.  We  have  found,  however, 
that  letting  letting  Rft  -  1 )  =  R[}(  —  1 )  =  a  I  eliminates  most  numerical  problems. 

One  Step  Ahead  Prediction 

I  he  n!h  order  least-squares  prediction  of  y(l)  ,  based  on  measurements  up  to  and  including  time  t-1, 
is 

n 

y„(0=£  y(‘~j)Anj(t-  1)  .  n  =  1 ,2 .  (3.12) 

y=i 


liquivalently,  the  least-squares  prediction  ofy’ft)  is  the  matnx  yn(t)  that,  if  taken  as  the  measurement  at 
time  t,  yields  e„(/)  =  0.  This  second  interpretation  of  ?„(/)  leads  to  an  algorithm  for  its  computation. 
Writing  (2.15)  in  the  notation  of  (2.25)  and  (2.54)  yields 


[/,(')]  =  [/-  />„«)]  t «A<0] 


()(mxp) 
«/'(•'-  1) 


(3.13) 


=  [/  -  O')]  [</>]>•(')  +  [/  -  O')] 


J 


The  fact  that  the  top  m  rows  of  [ fn{t)]  are  O')  and  30  argument  similar  to  the  argument  involving 
(2.47)-(2.50)  show  that  the  top  m  rows  of  (3.13)  are 


i 
i 
> 
, 1 

i 


( 


00=0'-  l)jt(')  +  O')  . 


where  en(t)  is  the  value  of  en{ i)  generated  by  the  lattice  for y(t)=  0.  The  prediction  then  is 


;■„(')=-  1)00  ,  n=  1.2  .... 


(3.14) 


This  prediction  can  be  generated  with  only  (3.2)-(3.4)  of  the  residual  error  lattice  algonthm. 


Including  Known  Inputs 

When  there  is  only  one  measurement  in  each  channel  (i.e.,  m=  1),  the  algonthm  presented  in  Section 
2  reduces  to  the  lattices  derived  in  [2]  and  [3].  One  of  the  major  advantages  of  the  algorithm  in  this 
paper,  as  mentioned  before,  is  its  capacity  to  accomodate  an  ARMA  model  m  which  all  outputs  have  the 
same  Auto  Regressive  structure.  For  such  a  model,  the  following  imbedding  accomodates  an  input  that 
does  not  have  the  same  MA  coefficients  for  each  channel. 


Let  there  be  m  scalar  outputs,  y}(i)  ,  and  one  scalar  input  u(t).  The  first  channel  then  is  [  rjU) . 

y}r,(t)  ]T  .  The  second  channel  is  the  m-vector  [  u(t),  0,...  0  ]'.  The  third  channel  is  [  U,  u(t),  0,...  0  ]', 
and  finally,  the  (m  +  1)  channel  is  the  m-vector  [  0,  0,  ...0,  u(t)  ] r 


It  is  easy  to  see  that  with  these  definitions,  the  desired  model  can  be  imbedded  into  a  p-channcl 
(p=m+  1)  vector-channel  Auto  Regressive  model.  The  extension  to  the  case  of  more  than  one  input  is 
obvious.  For  m  outputs  and  k  inputs,  the  resulting  vector-channel  lattice  has  one  channel  for  the  outputs 
and  mk  channels  for  the  inputs  (and  the  corresponding  copies). 


4  EXAMPLE 


Description  of  the  Structure  and  Data  Generation 


A  flexible  beam  is  cantilivered  to  a  rigid  disc  with  fixed  center,  and  a  torsional  spring  is  attached  to 


the  disc.  In-plane  motion,  including  the  linear  transverse  vibration  of  the  beam  is  modeled.  The  beam  is 
modeled  as  an  Euler- Bernoulli  beam  with  the  properties  listed  in  Table  1.  The  system  has  one  input,  the 
torque  applied  to  the  hub,  and  two  outputs,  the  elastic  tip  deflection  rj  and  hub  rotation  8. 


EIGURE  1.  The  Beam-Hub  Model 


TABLE  1.  Structural  Characteristics 


Mass  per  unit  length 
Hub  moment  of  inertia 
El  for  the  beam 
Torsional  spring  constant 
Beam  length 
Hub  radius 

Modal  Damping  Ratio 


=  .000288  slug/  inch 
=  85.6  slg  in2 
=  40,506,000  slug  in3/sec2 
=  85,600  lb-in 
=  1036.7  inches 
=  46  inches 

=  .005  (.5  %  critical)  in  each  mode 


To  generate  input-output  data,  the  state  space  representation  of  the  first  12  modes  of  the  structure 
was  excited  by  a  sequence  of  piecewise  constant  torques.  A  finite  element  model  with  15  uniform  beam 
elements  (31  degrees  of  freedom)  was  used  to  compute  the  first  12  modes  of  the  structure.  The  length  of 
the  sampling  interval  and  the  duration  of  the  piecewise  constant  torques  was  .05  seconds  (20  samples  per 
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second).  Hie  input  torque  sequence  was  obtained  from  an  IMSL  subroutine  that  generates  Gaussian 
random  sequences  with  mean  zero  and  variance  10,000. 


Fable  2  contains  the  true  discrete-time  system  eigenvalues  and  their  magnitudes,  and  the  damping 
ratios  (/eta)  and  frequencies  (Frcq  in  rad/s)  for  the  corresponding  continuous-time  eigenvalues.  (From 
here  on,  true  system  means  the  first  12  modes  of  the  finite  element  model.)  The  continuous-time 
frequencies  and  damping  ratios  were  obtained  from  the  logarithms  of  the  discrete-time  eigenvalues,  l  he 
first  five  modes  of  the  system  are  below  the  Nyquist  frequency  (10  Hz).  For  these  modes,  the  frequencies 
and  damping  ratios  in  the  first  five  rows  of  Fable  2  are  correct.  For  the  modes  above  the  Nyquist 
frequency,  aliasing  produces  frequencies  lower  than  the  correct  values  and  damping  ratios  higher  than  the 
correct  values.  (  The  aliased  frequencies  and  damping  ratios  arc  listed  in  Fable  2  for  comparison  with  the 
subsequent  eigenvalue  estimates.)  in  Table  2,  eigenvalues  are  ordered  according  to  magnitudes. 


TABLE  2.  True  System  Eigenvalues  (Discrete-  l  ime) 


Real  Fart  Imaginary  Fart 

/eta 

Freq 

Magnitude 

1  0.999390 

0.034733 

0.005000 

0.69 4825 

0.999827 

2  0.961521 

0.269716 

0.005000 

5.469626 

0.998634 

3  0.676489 

0.730860 

0.005000 

16.480264 

0.995888 

4  -0.099991 

0.986622 

0.005000 

33.435967 

0.991676 

5  -0.932683 

0.319972 

0.005000 

56.222164 

0.986043 

6  -0.444545 

0.872270 

0.010383 

40.842544 

0.979018 

7  0.924166 

0.296545 

0.095783 

6.209997 

0.970578 

8  -0.164358 

0.946500 

0.023021 

34.854594 

0.960664 

9  -0.515252 

0.797210 

0.024291 

42.891594 

0.949225 

10  0.766445 

0.537675 

0.107091 

12.234877 

0.936233 

11  -0.759660 

0.521941 

0.032095 

50.791924 

0.921686 

12  0.503091 

0.753000 

0.100486 

19.635824 

0.905600 

For  studying  the  effects  of  measurement  noise,  small  amounts  of  white  noise  were  added  to  the 
outputs  and  the  resulting  noisy  measurements  were  fed  into  the  lattice  filters.  Note  that  three  uncorrelated 
sequences  of  white  noise  were  used.  The  first  was  the  torque,  which  was  an  input  known  to  the  lattice; 
hence  this  sequence  was  white  noise  only  in  the  sense  of  its  statistics  and  the  lattice  saw  it  as  a  known 
deterministic  input.  The  other  two  white  noise  sequences  were  scaled  by  various  amounts  and,  in  some 
cases,  added  to  the  data  from  the  two  sensors  to  produce  the  noisy  measurement  sequences  that  were  fed 
into  the  lattice.  This  measurement  noise  was  not  known  to  the  lattice. 


It  is  well  known  that,  in  the  presence  of  measurement  noise,  the  least-squares  method  yields  biased 
estimates  for  the  AR  coefficients.  If  it  is  essential  to  obt;un  accurate  estimates  for  these  parameters,  other 
identification  methods  such  as  extended  least  squares  (ELS)  or  instrumental  variable  (IV)  methods  can  be 
used.  The  least-squares  lattice  can  be  modified  easily  for  the  ELS  method.  For  lattice  implementation  of 
the  IV  method,  see  [15]  or  [16]. 

We  processed  input/output  data  with  the  lattice  filters  to  estimate  the  number  of  excited  modes  (one 
half  of  the  plant  order)  and  the  ARMA  coefficients.  From  the  AR  coefficients,  we  obtained  estimates  for 
plant  eigenvalues.  We  used  a  vector-channel  lattice  and  two  standard  lattices  for  comparison.  The 
vector-channel  lattice  had  three  channels,  each  with  two  components,  (p  =  3,  m  =  2).  Channel  1  contained 
the  two  measurements  (hub  rotation  and  tip  deflection);  Channels  2  and  3  contained,  as  described  at  the 
end  of  section  3,  the  required  copies  of  the  input  torque.  We  used  the  two  standard  lattices  independently 
for  the  hub  rotation  measurement  and  for  the  tip  deflection  measurement.  Each  of  these  lattices  had  two 
one  dimensional  channels  (p=2,  m=  1),  Channel  1  for  the  measurement  (hub  or  tip)  and  Channel  2  for 
the  input  torque.  Due  to  the  ambiguous  relationship  between  minimal  lattice  order  and  the  number  of 
excited  modes,  mentioned  in  Section  i,  and  the  possibility  of  extraneous  eigenvalues,  we  did  not  use  both 
measurements  together  in  a  standard  multi-input-multi-output  lattice.  For  all  runs,  the  value  of  the 
foreettinu  factor  /.  was  0.95. 


Order  Determination 


Recall  the  definition  (2.23)  for  R'n(t)  .  Each  diagonal  element  of  this  matrix  is  the  square  of  the  n!h 
order  forward  error  norm  for  the  corresponding  channel.  Figure  2,  which  was  obtained  from  the 
vector-channel  lattice  after  300  steps,  shows  the  behavior  of  the  (1,1)  clement  of  A’‘'(f)  (the  norm  of  the 
forwasr  error  for  the  first  vector  measurement  channel)  as  n  increases.  Recall  that  the  first  channel  of  the 
vector-channel  lattice  contains  the  measurements  from  both  sensors.  The  lowest  curve  corresponds  to  the 
no-noisc  case.  The  other  curves  correspond  to  successively  increasing  amounts  of  noise  added  to  the 
measurements. 

Ibe  effective  order  of  the  structure,  which  is  twice  the  number  of  significantly  excited  modes,  can  be 
determined  from  how  R%(t)  varies  with  the  lattice  order  n  The  most  important  feature  of  the  three  curses 
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in  Figure  2  is  the  drop  as  the  lattice  order  goes  from  10  to  12  (i.e..  from  5  to  6  modes).  1'his  drop 
indicates  that  the  effective  order  of  structure  is  12. 

The  behavior  of  the  (1,1)  element  of  A’,'j(/)  in  the  absence  of  measurement  noise  indicates  that  six 
modes  are  excited  substantially  by  the  input.  The  remaining  six  modes  in  the  simulation  model  are 
excited  minimally  and  contribute  little  to  the  output.  The  drastic  drop  in  the  graph  between  n=  10  and 
n  =  12  occurs  when  the  order  of  the  fdter  matches  the  number  ot  substantially  excited  modes.  (  Ihe  plant 
order  must  be  even.)  As  n  increases  past  12,  more  and  more  of  the  output  from  the  minimally  excited 
modes  is  accounted  for,  until  the  forward  error  falls  to  /.ero  at  n  -  24  because  a  24'/l  order  ARM  A  model 
tits  all  of  the  output. 

When  small  amounts  of  noise  are  added  to  the  measurements,  the  graph  changes.  I  he  drop  still 
occurs  between  10  and  12,  but  its  magnitude  decreases  as  the  amount  of  measurement  noise  increases. 
Furthermore,  as  n  increases  past  12,  the  (1,1)  element  of  A IJjU)  hardly  decreases  at  all.  The  excitation  of 
the  structure  is  the  same  with  and  without  measurement  noise,  but  the  output  from  the  minimally  excited 
modes  is  lost  in  the  measurement  noise  and  the  filter  cannot  distinguish  these  modes  from  the  noise. 
With  increasing  measurement  noise,  the  graph  approaches  a  level  straight  line,  indicating  that  the  filter  is 
distinguisliing  less  and  less  data  from  the  structure  and  noise. 

The  accuracy  of  the  subsequent  eigenvalue  estimates  indicates,  as  expected,  that  the  five  sub-Nyquist 
frequency  modes  are  excited  most  and  the  first  mode  above  the  Nsquist  frequency  is  excited  substantially 
though  lass  than  the  first  five.  Different  input  sequences,  including  combinations  of  many  sine  waves, 
produced  the  same  behavior.  In  all  cases,  only  the  frequencies  below  Nyquist  frequency  are  excited 
substantially.  This  is  not  surprizing,  since  the  input  is  applied  at  the  sampling  frequency. 

.Several  authors  have  suggested  criteria  to  estimate  the  order  of  a  system  from  input-output  data,  see 
[9,17],  Fach  of  these  criteria  looks  for  a  minimum  of  a  scalar  function  of  the  form 

F(n)=i  in  [  1  It  y  K(0]  +  " */(')•  (4.1) 

The  major  difference  among  these  criteria  is  the  choice  of/(t)  .  For  the  systems  studied  here,  a  variety  of 
these  methods  were  attempted  and  none  performed  noticeably  better  than  the  (1,1)  clement  of  A£(f)  alone, 
as  used  in  Figure  2. 
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To  estimate  tlie  eigenvalues  of  the  structure,  wc  generated  the  estimates  of  the  AR  coefficients  with 
the  algorithm  in  Section  3,  for  various  orders,  and  computed  the  eigenvalues  of  the  AR  model.  The  only 
case  in  which  the  estimated  AR  coefficients  were  correct  was  where  there  was  no  sensor  noise  and  the 
lattice  order  was  24  (full  order;  i.e.,  the  order  of  the  true  system),  but  the  estimated  eigenvalues  for  the 
sub-Nyquist  modes  are  generally  good  for  low  measurement  noise  and  lattice  order  18  or  greater.  This  is 
possible  because  most  of  the  AR  coefficients  of  a  large  order  system  are  highly  sensitive  to  perturbations 
in  the  eigenvalues  and,  conversely,  the  eigenvalues  arc  relatively  insensitive  to  perturbations  in  the  AR 
coefficients.  The  relationship  established  by  Cayley- Hamilton  between  the  AR  coefficients  and  the 
eigenvalues  makes  the  high  sensitivity  of  AR  coefficients  with  respect  to  eigenvalues  apparent. 

In  tins  section,  we  are  interested  primarily  in  the  effects  of  measurement  noise  and  unmodeled  inodes. 
In  the  absence  of  noise,  till  the  eigenvalues  corresponding  to  the  higher  modes  were  identified,  given  high 
enough  filter  order.  The  full-order  vector-channel  lattice  (order  24)  identified  all  of  the  system  eigenvalues 
w  ith  four  or  more  significant  digits,  except  one  of  the  faster  modes,  which  had  only  three  digit  accuracy 
The  full-order  standard  lattice  for  either  measurement  produced  eigenvalue  estimates  with  2  to  3  accurate 
significant  digits.  Table  3  compares  standard  and  vector-channel  lattices  for  the  18,A -order  filter  with  no 
noise  added  to  the  measurements. 

In  Tables  3  and  4,  complex  eigenvalues  arc  represented  by  their  continuous-time  damping  ratio  and 
trequency  and  discrete-time  magnitude.  The  frequencies,  damping  ratios  and  magnitudes  estimated  for  the 
tirst  five  modes  are  listed  first.  Then,  the  aliased  frequencies  and  damping  ratios  estimated  for  the 
remaining  modes  are  ordered  according  to  the  magnitude  of  the  corresponding  discrete-time  eigenvalues, 
fil  hesc  abased  frequencies  and  damping  ratios  are  all  that  can  be  identified  from  the  input/output  data, 
without  other  knowledge  of  the  true  system.)  In  practice,  the  magnitudes  of  the  eigenvalues  arc  bkely  to 
be  the  only  tool  for  proper  ordering.  [t  is  common  in  structures,  for  example,  that  discrete-time 
eigenvalues  corrsponding  to  higher  modes  have  smaller  magnitudes. 

To  study  the  effects  of  measurement  noise,  v  ■  added  small  amounts  of  noise,  roughly  3%  of  the 
output  on  average,  to  the  measurements.  Table  4  compares  the  scalar  and  vector-channel  lattices  for  this 


TABLE  3.  lattice  Comparisons  :  No  Noise 

Order  = 

18 

After  1500  Steps 

Zeta 

Ereq 

Mag 

l  ip  Disr 

dacement  Lattice 

-0.1793 

0.6866 

1.0063 

0.0057 

5.4705 

0.9984 

0.0050 

16.4786 

0.9959 

0.0051 

33.4305 

0.9914 

0.0050 

56.2204 

0.9861 

0.0122 

40.7381 

0.9754 

0.0474 

48.4888 

0.8913 

0.0785 

34.8581 

0.8717 

positive  real 

0.2882 

negative 

real 

0.1907 

Hub  Rotation  Lattice 

0.1301 

0.6933 

0.9955 

0.0060 

5.4721 

0.9984 

0.0050 

16.4791 

0.9959 

0.0052 

33.4310 

0.9913 

0.0050 

56.2205 

0.9861 

0.0128 

40.7135 

0  9744 

0.0487 

48.6059 

0.8882 

0.0953 

35.2980 

0.8445 

0  8384 

31.6749 

0.0875 

Vector-Channel  lattice 

0.0809 

0.6987 

0.9972 

0.0060 

5.4696 

0.9984 

0.0050 

16.4790 

0.9959 

0.0052 

33.4304 

0.9914 

0.0050 

56.2214 

0.9861 

0.0125 

40.7320 

0.9748 

0.0470 

48.6657 

0.8919 

0.0886 

35.1743 

0.8551 

positine  real 

0.1289 

negative  real 

0.0697 
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TABLE  4.  Lattice  Comparison:  Noisy  Measurements 


Order  =  24  (Full  Order) 
After  1200  Steps 

Zeta  Freq  Mag 
Tip  Displacement  Lattice 


Order  =  20 
After  1 500  Steps 
Zeta  Freq  M 

Tip  Displacement  Lattice 


Order  =  1 8 
After  1 500  Steps 
Zeta  Freq  5 

Tip  Displacement  lattice 


-0.0813 

0.7219 

1.0029 

0.9710 

0.2177 

0.9567 

0.0031  5.4591 

0.0046 

5.4668 

0.9988 

0.0031 

5.4718 

0.9992 

0.0040  16.4651 

0.0052 

16.4814 

0.9957 

0.0040 

16.4631 

0.9967 

0.0057  33  4285 

0.0056 

33.4273 

0.9907 

0.0057 

33.4298 

0.9906 

0.0041  56.3500 

0.0056 

56.2121 

0.9845 

0.0039 

56.3435 

0.9890 

-0.0059  23.8576 

0.0103 

40.7798 

0.9793 

0.0141 

23.6348 

0.9835 

0.0133  40  8553 

0.0129 

44.9677 

0.9714 

0.0119 

40.8653 

0.9760 

0.0203  46.4778 

0.0151 

50.1422 

0.9627 

0.0142 

46.7545 

0.9674 

positive  real 

0.0676 

22.8835 

0.9254 

0.4058 

24.7266 

0.5776 

negative  real 

0.0564 

28.5524 

0.9226 

negative  real 

0.9338 

positive  real 

0.3655 

15.4049 

0.7390 

negative  real 

0.3413 

positive  real 

negative  real 
negative  real 

Hub  Rotation  Lattice 


0.9540 

0.3509 


16.4795 

33.4165 

56.2169 

40.7037 

45.4384 

51.3041 

27.9494 

25.0789 


negative  real 
negative  real 


Hub  Rotation  Lattice 


Hub  Rotation  lattice 


0.0593 

0.6802 

0.7049 

0.0153  5  4965 

0.99: 

>8 

0.0048 

5.4703 

0.9987 

0.0108 

5.4909 

0.9970 

0  0044  16  4938 

0  9964 

0.0051 

16.4789 

0.9958 

0.0049 

16.4915 

0.9960 

0.0054  334298 

0.9909 

0.0054 

33.4057 

0.9911 

0.0054 

33.4236 

0.9909 

0.0064  56.2400 

0  9821 

0.0052 

56.2204 

0.9854 

0.0059 

56.2616 

0.9835 

0.0014  41  1753 

0.9972 

0.0035 

45.9140 

0.9921 

-0.0057 

48.1314 

1.0138 

0.0527  24  0119 

0.9386 

0.0123 

40.6022 

0.9754 

0.0048 

41.0987 

0  9901 

0.0282  486725 

0.9335 

0.0192 

51.8548 

0.9515 

0.0951 

21.7632 

0.9013 

positive  real 

0  9648 

0.0668 

27.8271 

0.9110 

0  1231 

27.0668 

0.8455 

positive  real 

0  9510 

0.0717 

26.0990 

09104 

0  0909 

58.5534 

0.7655 

negative  real 

0.7082 

0.1526 

15.2101 

0.8892 

negative  real 

0  1385 

negative  real 

0.9563 

negative  real 

0.2579 

Vector-Channel  Lattice 

Vector-Channel  Lattice 

Vector-Channel  Lattice 

0.0133 

0.6941 

0.9995 

0.6853 

1 

0.0044 

5.4658 

0.9988 

0.0052 

5.4792 

0.9986 

U.99i 

0.0044  It 
0.0055  3: 

0.0055  5< 

0.0044  41 

0.0074  M 
0.0553  2; 

0.2774  2< 

negative  real 
negative  real 


16.4813 

33.4307 

56.2822 

48.2454 

40.9636 

23.5206 

26.0845 


0.9964 

0.9908 

0.9848 

0.9894 

0.9849 

0.9370 

0.6863 

0.8681 

0.4307 


0.0043 

16.4825 

0.0056 

33.4308 

0.0056 

56.2721 

0.0080 

40.9916 

0.0193 

24.1705 

0.0174 

48.0543 

negative 

real 

negative 

real 

>33 
i ;<a 


1 


Overall,  the  lattices  produced  good  estimates  for  the  modes  below  the  Nyquist  frequency.  The 
eigenvalue  corresponding  to  the  fust  mode  above  the  Nyquist  frequency  also  was  obtained  with 
reasonable  accuracy,  although  the  frequency  and  the  damping  ratio  for  this  mode  arc  aliased.  The  lowest 
mode,  corresponding  to  the  hub  rotation,  is  the  hardest  among  the  first  five  to  identify.  Longer  data  sets 
and  higher-order  filters  are  needed  to  obtain  reasonable  estimates  for  this  mode.  In  general,  the 
reduced-order  lattices  require  more  data  points  to  obtain  good  estimates  of  frequencies  and  damping 


Generally,  lattices  of  order  lower  than  18  missed  the  lowest  mode  entirely,  but  identified  the  other 
four  sub-Nyquist  frequencies  satisfatonly,  although  some  of  the  damping  ratios  w'ere  quite  inaccurate.  If 
we  only  need  small  prediction  error  from  the  filter  (say,  for  adaptive  control),  then  figure  2  indicates  that 
12  is  a  sufficient  lattice  order. 

When  we  added  sensor  noise  that  was  on  average  5%  of  the  sensor  data,  the  results  were  almost  as 
good  as  those  in  the  tables  for  3%  noise,  for  noise  levels  of  9%  and  above,  the  filter  identified  only  one 
or  two  of  the  higher  sub-Nyquist  frequencies  and  yielded  poor  estimates  for  damping  ratios.  This  is 
consistent  with  figure  2,  wliich  indicates  that  the  lattice  filter  has  trouble  distinguishing  between  the 
structural  dynamics  and  the  greater  sensor  noise.  We  expect  that  the  estimates  of  the  lowest  frequencies 
would  improve  with  more  samples,  which  would  cover  more  periods  of  the  lower  frequencies. 

I  ables  3  and  4  indicate  the  advantages  of  the  vector-channel  lattice  The  vector-channel  estimates  for 
the  lowest  mode  are  significantly  better  than  those  obtained  by  the  independent  standard  lattices,  for  the 
rest  ot  the  modes,  where  standard  lattices  produced  good  estimates,  the  vector-channel  produced  estimates 
close  to  the  average  of  the  two  independent  standard  lattices. 

The  ordenng  of  the  eigenvalues  is  also  important.  Often  the  standard  lattices  estimated  erroneous 
eigenvalues  whose  magnitudes  were  greater  than  the  magnitudes  of  the  true  eigenvalues.  Sometimes  these 
erroneous  eigenvalues  even  had  magnitudes  greater  than  one  (and  negative  damping  ratios).  This  greatly 
complicates  identifying  the  true  modes  below  the  Nyquist  frequency. 

We  found  that,  in  a  large  majority  of  comparisons,  the  vector-channel  lattice  estimates  were  at  least 
as  good  as  or  better  than  the  estimates  from  either  standard  lattice.  However,  since  the  vector-channel 
estimate  is  a  kind  of  average  of  the  individual  standard  lattice  estimates,  it  is  possible  that  one  sufficiently 
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noisy  measurement  channel  can  cause  the  vector-channel  estimate  to  be  worse  than  another  single  channel 
estimate.  litis  happened  in  Table  4.  for  example,  wheie  the  standard  lattice  for  the  hub  channel  estimated 
a  negative  damping  ratio,  which  caused  the  vector-channel  lattice  to  estimate  an  incorrect  mode  with 
frequency  48. 2  and  damping  ratio  similar  to  that  of  the  true  modes. 


s  CONTI  l  SIONS 


We  have  developed  a  least-squares  lattiee  filter  that  can  constrain  the  coefficients  in  the  AR  model  to 
be  the  same  lor  some  or  all  of  the  channels  of  a  multi-output  system.  If  this  constraint  is  not  imposed, 
the  vector-channel  lattice  here  reduces  to  the  standard  multi-input-multi-output  least-squares  lattice  filter, 
liven  in  the  standard  lattice  case,  the  derivation  in  Section  2  is  needed  for  current  research  on  the 
performance  of  lattice  filters  applied  to  distnbuted  systems,  which  require  infinite  histories  in  the  AR 
models,  l  o  handle  forced-response  data,  as  opposed  to  impulse  response,  the  vector-channel  lattice  must 
accomodate  scalar  input  (MA)  channels  as  well  as  vector  measurement  channels.  A  method  for  this  has 
been  demonstrated,  although  the  number  of  channels  grows  large  fast  w  ith  multiple  inputs. 

The  primary  motivation  for  the  vector-channel  lattice  has  been  adaptive  identification  and  control  of 
flexible  structures.  Constraining  all  measurement  channels  to  have  the  same  AR  model  introduces 
information  into  the  estimation  algorithm  that  reduces  the  effect  of  sensor  noise,  as  demonstrated  In  the 
numerical  results  in  Section  4.  The  numerical  results  also  show  that  the  frequencies  and  damping  ratios 
tor  the  most  significant  modes  can  be  identified  with  a  lattice  whose  order  is  lower  than  the  number  of 
excited  modes,  even  in  the  presence  of  noise. 
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Abstract 

In  this  paper,  we  present  our  recent  investigations  on 
lattice  filters  and  their  applications  to  adaptive  identifica¬ 
tion  of  flexible  structures.  Since  the  order  of  the  systems 
can  not  be  known  -or  the  effective  order  may  change-  the 
order  recursiveness  of  the  lattices  is  of  particular  inter¬ 
est.  Implementation  of  lattices  would  permit  on-line  or¬ 
der  determination  and  would  allow  the  order  of  the  filter 
to  be  changed  without  the  need  to  reprocess  the  previous 
data.  Experimental  data  from  the  flexible  grid  structure 
at  NASA-Langley  are  used  to  obtain  results  showing  the 
feasibility  of  lattices  and  the  advantages  that  result  from 
their  order  recursive  property.  One-step-ahead  prediction 
and  estimates  for  natural  frequencies  are  among  the  results 
shown.  Of  particular  interest  are  the  frequency  estimates 
which  agree  cloaely  with  the  frequency  estimates  obtained 
from  off-line  identification  techniques.  The  one  step-ahead- 
prediction  results  also  show  the  advantages  that  lattices 
provide  with  their  order-determination  capability,  which 
would  be  significant  for  adaptive  control  purposes. 


Montgomery,  Sundararajan  [in  SMl,SM2]  and  Wiberg 
[Wl]  were  first  to  apply  lattices  to  identification  of  flexible 
structures.  Montgomery,  et.  al.,  have  been  using  a  variety 
of  techniques  -such  as  Discrete  Fourier  Tranforms-  for  the 
identication  of  these  structures  while  the  lattices  were  re¬ 
lied  upon  mostly  for  order  determination.  Wiberg  has  been 
emphasizing  mainly  the  frequency  estimation  based  on  the 
free  response  (zero  input)  of  the  system.  In  our  research 
|JG1,J1],  we  have  expanded  the  scope  and  developed  new 
approaches  in  the  application  of  lattice  filters  for  adaptive 
identification  of  flexible  systems.  Our  approach  uses  the 
lattice  filters  as  the  only  identification  tool  required.  We 
also  seek  frequency  (and  damping  ratio)  estimates  for  a  va¬ 
riety  of  input  sequences  with  the  eventual  goal  of  adaptive 
control. 

In  this  paper,  we  present  results  obtained  from  the  ap¬ 
plication  of  lattice  filters  to  the  NASA-Langley  grid,  a  flexi¬ 
ble  experimental  structure  used  in  research  on  control  and 
identification  of  space  structures.  These  results  demon¬ 
strate  the  feasibility  of  lattice  filters  for  on-line  identifica¬ 
tion  -  and  eventually  adaptive  control  -  of  flexible  struc- 


Adaptive  identification  of  flexible  systems  is  receiving 
increasing  attention.  Since  flexible  systems  are  often  mod¬ 
eled  as  distributed  systems,  they  provide  major  practical 
and  theoretical  challenges  to  the  control  community.  Insuf¬ 
ficient  model  fidelity,  poor  knowledge  of  structural  damp¬ 
ing,  and  the  lack  of  feasibility  for  ground  testing  for  many 
of  the  large  space  structures  are  some  of  the  contributing 
factors  to  the  recent  interest  in  adaptive  identification  of 
these  systems.  Determination  of  an  effective  system  order 
is  one  of  the  most  important  problems,  since  there  is  no 
fixed  order  that  can  be  assigned  to  these  systems  a  priori. 

Recently,  lattice  filters  are  being  studied  for  adaptive 
identification  of  flexible  structures  because  they  can  iden¬ 
tify  the  effective  plant  order  as  well  as  the  parameters  for 
the  digital  input/output  model  of  a  structure.  The  time- 
update  equations  needed  for  adaptive  filtering  were  de¬ 
veloped  by  Lee,  Morf,  and  FYeidlander  who  studied  the 
applications  of  lattices  to  speech  processing  and  adaptive 
control,  see  [LMFJLFM]  among  many. 

‘This  uUor  was  •apportad  by  AFOSR  put  (40109 


2.  State-space  and  ARMA  representations 

For  numerical  work,  the  linear  motions  of  flexible  struc¬ 
tures  are  approximated  by  a  finite  dimensional  differential 
equation  of  the  form 

f  i(t)  =  Ai(t)  +  Bcu(t) 

1  V(t)  -  Cx(t ) 

where  y  is  the  output  (measurements),  u  the  input,  etc. 
The  discrete-time  form  of  the  above  equation  has  the  fol¬ 
lowing  well  known  representation: 


f  x(t  +  1)  =  Tx(t)  +  B4u(t) 
\  y(t)  =  Cx{t) 


For  digital  adaptive  control  and  identification,  the  input 
-  output  representation  corresponding  to  (2.2)  usually  is 
used.  This  alternative  representation  is 


v(f)  =  Z!  «.  v(t  - 1)  +  £  «(*  -  0  (2-3) 

im  I  i«I 
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and  is  called  an  ARMA  (Auto  Regressive  Moving  Average) 
model. 

Mass  and  stiffness  properties  of  a  system  are,  tradition¬ 
ally  speaking,  obtained  from  finite  element  analysis  and  ex¬ 
tensive  ground  testing.  The  matrices  in  equations  (2.1,  2.2) 
could  then  be  formed  by  the  standard  methods.  This  ap¬ 
proach,  as  mentioned  in  the  previous  section,  is  becoming 
less  and  less  desirable,  due  to  the  properties  and  require¬ 
ments  associated  with  the  modem  applications  of  flexible 
systems. 

Because  of  the  difficulties  associated  with  the  identifi¬ 
cation  of  the  matrices  in  (2.1)  and  (2.2)  -  such  as  the  num¬ 
ber  of  parameters  involved  and  the  lack  of  uniqueness  -  the 
ARMA  representation  of  (2.3)  is  more  useful  for  adaptive 
identification  and  control.  In  addition  to  the  small  number 
of  parameters  involved,  this  model  uses  only  the  measured 
inputs  and  outputs  of  the  system  (the  state  need  not  be 
estimated)  and  it  can  be  used  for  a  variety  of  adaptive 
control  algorithms  (see  [GSlj  among  many).  Furthermore, 
the  AR  coefficients  contain  the  necessary  information  for 
natural  frequency  (and  damping)  determination.  This  can 
be  seen  easily  by  considering  the  single-input-single-  out¬ 
put  (SISO)  case.  There,  the  AR  coefficients  are  the  same 
as  the  coefficients  of  the  characteristic  polynomial  of  the 
state  transition  matrix  of  equation  (2.2)  (for  details  see 
[JG1  and  Clj).  In  [JG1,  J 1]  we  have  presented  methods 
that  extends  this  directly  to  the  MIMO  case,  as  well. 

Adaptive  (on-line)  identification  of  these  ARMA  mod¬ 
els  has  been  studied  extensively  by  many  authors,  see 
[LSI,  GSlj  for  example.  In  most  cases,  however,  the  order 
of  the  system  is  assumed  to  be  known  and  fixed.  Such 
assumptions,  clearly,  are  not  applicable  to  the  flexible  sys¬ 
tems  of  interest.  Lattice  filters,  introduced  in  the  next 
section,  have  the  capability  to  identify  the  required  param¬ 
eters  without  such  assumptions  and  can  actually  estimate 
the  effective  order  of  the  system  as  well. 

3.  Adaptive  Identification  Using  Lattice  Filters 

The  most  common  method  for  the  on-line  estimation 
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Table  1:  The  Residual  Error  Filter  Algorithm 

all  intermediate  order  filters  are  automatically  obtained, 
the  effective  order  of  the  system  can  be  estimated  so  that 
the  most  appropriate  set  of  coefficients  (corresponding  to 
the  best  estimated  order)  are  used. 

The  order  determination  capability  of  lattices  is  the 
major  reason  for  our  interest.  It  should  be  noted,  how¬ 
ever,  that  other  important  features  of  lattices  are  their 
numerical  stability  [LLl]  and  their  potential  for  VLSI  im¬ 
plementation  due  to  their  special  structure  which  consists 
of  many  similar  modules. 

In  [JGl,Jl),  we  have  developed  a  novel  derivation  of 
the  lattices  that  has  two  important  features.  The  first  one 
is  the  basic  framework  used,  which  lends  itself  nicely  to 


of  the  ARMA  coefficients  is  the  recursive  least-squares  al¬ 
gorithm.  In  most  applications,  the  fact  that  N  is  finite 
is  not  a  serious  limitation  because  the  dynamic  behavior 
of  the  system  can  be  approximated  by  the  behavior  of  a 
finite  dimensional  approximation  to  the  system.  For  flex¬ 
ible  structures,  which  theoretically  speaking  are  infinite 
dimensional,  it  is  quite  common  to  use  an  approximation 
model  based  on  the  few  of  the  modes  that  are  highly  ex¬ 
cited.  The  real  problem  is  that,  in  the  classical  recursive 
least-squares  algorithms,  N  must  be  fixed.  Depending  on 
the  initial  conditions  and  external  excitation,  the  order  of 
the  finite  dimensional  model  required  to  approximate  the 
response  of  the  system  accurately  can  vary  widely. 

The  least-squares  lattice  addresses  this  issue  directly 
by  providing  an  algorithm  that  can  be  variable-order.  A 
lattice  filter  algorithm  of  order  n,  being  recursive  in  time 
as  well  as  order,  results  in  filters  of  order  one  to  n  and 
allows  the  maximum  order  of  the  filter  (n)  to  be  increased 
gradually,  without  reprocessing  of  the  previous  data.  Since 


some  of  the  more  theoretical  aspects  of  these  flexible  sys¬ 
tems.  The  second  feature  is  the  extension  of  the  lattices  to 
the  case  of  many  measurement  channels  (where  each  chan¬ 
nel  may  contain  many  scalar  measurements).  This  feature 
can  be  exploited  easily  to  extend  the  many  desirable  prop¬ 
erties  of  the  SISO  ARMA  models  to  the  MIMO  systems 
and  hence  improve  the  estimates  for  the  system  natural 
frequencies,  see  [JGl,Jlj  for  more  detail. 

The  residual-error  lattice  algorithm  is  listed  in  Table 
I.  This  algorithm  must  be  performed  at  each  time  step  to 
update  the  filter  variables.  The  variable  R*(t)  is  one  of  the 
crucial  variables  and  is  used  to  help  determine  the  effective 
order  of  the  system.  Due  to  space  limitation,  the  discussion 
concerning  this  issue  is  omitted.  For  details  regarding  the 
initialization  and  definitions  of  other  variables  see  [JGl,JlJ. 

The  algorithm  in  Table  1,  as  mentioned  before,  calcu¬ 
lates  the  residual  error;  i.e.,  the  portion  of  the  output  y(t) 
that  can  not  be  fitted  in  a  model  based  on  all  of  the  data 
including  y(t).  A  more  useful  approach  (for  prediction  and 
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Table  2:  AR  Coefficients  Algorithm 

eventually  control  purposes)  would  be  the  prediction  error 
filter.  By  prediction  error  we  mean  the  error  between  y(t) 
and  its  estimate  using  the  data  up  to  time  (t-l).  Simple 
modification  of  the  above  algorithm  would  yield  a  filter 
that  calculates  the  prediction  error  (and  eliminates  the 
inversion  of  the  G  matrix  in  the  process).  Also,  simple 
application  of  the  matrix  inversion  formulas  results  in  the 
elimination  of  another  inversion  (preferably  R')  and  thus 
reduce  the  total  number  of  inversions  to  one  for  each  time 
step.  Again,  we  omit  the  details  for  the  sake  of  brevity. 
The  prediction  error  lattice  algorithm  would  be  quite  sim¬ 
ilar  to  the  residual  algorithm.  Four  of  the  above  equations 
have  to  be  modified  as  the  following  (  with  e  and  f  as  the 
prediction  errors  )  : 

Kn+l(t)  =  X  K„+i(t  -  1)  +  e„(t)  G„(f  -  1)  r„(t  -  1) 
Gn+1(t)  =  G„(<)  +  G„(t)  f„(t)  R;’(t  -  1)  f„(t)  G„(t) 

Mi(0  =  MO  -  M«  -  i)  R;r(‘  -  2)  /C,(‘  -  1) 
M,(0  =  M<  -  1)  -  MO  K'(t  -  0  *.♦,(*  -  1) 

Table  2  consists  of  the  algorithm  that  generates  the 
ARMA  coefficients  for  any  desired  order.  For  this,  the 
ARM  A  of  (2.3)  is  first  converted  to  a  two-channel  AR  form 
(using  a  simple  embedding).  This  part  is  presented  sepa¬ 
rately  to  underscore  the  point  that  these  coefficients  need 
not  be  updated  at  every  time  step.  Rather,  this  algorithm 
is  invoked  at  any  time  (and  for  any  order)  desired  to  gener¬ 
ate  the  ARMA  coefficients  and  will  be  turned  off  at  other 
times. 

Other  technical  issues  (such  as  the  required  embedding, 
the  reduced  lattice  structure  to  determine  the  input  for  the 
control  applications,  etc.)  are  discussed  in  (Jlj. 
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Figure  1:  The  Grid  Schematic 


4.  The  NASA-Langley  Grid 

This  grid  is  an  experimental  structure  at  the  NASA- 
Langley  Research  Center  in  Hampton,  Virginia,  designed 
for  research  on  identification  and  control  of  flexible  space 
structures.  It  is  a  7  ft  by  10  ft  planar  structure  made  of 
thin,  flat  aluminum  bars  with  uniform  spacing  (8  vertical 
and  11  horizontal  bars).  The  grid  is  suspended  from  the 
ceiling  by  two  cables  as  shown  in  the  schematic  of  Figure  1. 
Rate  gyros  are  used  to  measure  the  rotational  rate  about 
an  axis  in  the  plane  of  the  grid  and  the  actuators  are  reac¬ 
tion  wheels.  The  position  of  the  actuators  and  sensors  are 
marked  on  the  diagram  below. 

This  structure  has  been  studied  extensively  by  Mont¬ 
gomery,  et  ai,  [MWLN.SMl].  Sine  wave  sweeps  have  been 
used  to  determine  the  structural  natural  frequencies  and 
the  results  show  that  the  grid  has  as  many  as  twenty  modes 
below  the  frequency  of  10  Hz.  These  closely  packed  natural 
frequencies  (and  their  corresponding  low  damping  ratios) 
are  the  main  reasons  for  our  interest  in  this  structure  for 
our  research  in  adaptive  identification  and  control  of  flex¬ 
ible  systems. 

5.  Experimental  Results 

In  this  section,  we  present  some  of  the  results  obtained 
from  the  grid  by  the  lattice  filter  algorithms  listed  above. 
All  of  he  results  presented  here  are  obtained  from  two  test 
runs.  In  the  first  case  the  input  consists  of  a  sine-wave 
with  a  varying  frequency.  The  second  case  has  a  square- 
wave  as  the  input  and  is  used  to  excite  numerous  modes 
of  the  structure.  In  each  case,  the  aim  is  to  identify  the 
most  important  natural  frequencies  (and  their  correspond¬ 
ing  damping)  and  study  the  behavior  of  the  prediction  er¬ 
ror  for  different  filter  orders. 

Figure  2  is  a  sample  of  the  sensor  measurements  for 
the  case  of  square  input.  The  top  plot  is  the  first  sensor, 
the  second  plot  corresponds  to  sensor  number  2,  etc.  It 
is  clear  that,  although  all  sensors  have  certain  amount  of 
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Figure  2:  Sample  of  Sensor  Measurements:  Sensors  1-4 


noise,  sensor  four  is  quite  noisy  and  sensor  three  seems  to 
be  saturated  most  of  the  time.  Overall,  sensor  two  has  the 
best  measurements. 

The  first  set  of  results  concern  the  prediction  error  and 
its  behavior  as  the  order  of  the  filter  increases.  Since  most 
of  the  available  adaptive  control  algorithm  are  based  on 
some  prediction  for  the  future  outputs  (one  or  more  steps 
ahead),  this  study  is  the  first  step  toward  the  eventual 
adaptive  control  of  these  systems.  Figue  4  shows  the  ac¬ 
tual  outputs  (the  first  graph)  and  the  (one  step  ahead)  pre¬ 
dicted  value  of  the  output  for  different  order  filters.  The 
saw-tooth  curve  around  the  horizontal  axis  is  the  actual 
prediction  error  for  that  order  filter.  The  error  is  plot¬ 
ted  on  the  same  scale  so  that  it  can  be  compared  with 
the  actual  prediction.  Also  note  the  drop  in  the  average 
magnitude  of  the  error,  as  the  filter  order  is  increased. 

As  these  graphs  show,  the  general  shape  of  the  response 
(correponding  to  the  dominant  modes)  is  identified  and 
predicted  by  relatively  low  order  lattices.  The  high  fre¬ 
quency  content  of  the  signal,  however,  requires  filters  of 
significantly  larger  order.  This  is  evident  in  the  peaks  of 
the  signal  where  the  high  frequencies  are  more  distinguish¬ 
able.  The  next  graph  (Figue  5)  is  the  magnified  version  of 
Figure  4.  There,  the  behavior  of  the  prediction  (and  the 
corresponding  error)  is  underscored  more  clearly. 

Next,  we  will  summarize  some  of  results  regarding  the 
natural  frequency  estimates  of  the  system.  We  have  cho¬ 
sen  to  present  the  results  from  the  case  of  sinewave  input. 
The  results  obtained  from  this  case  show  that  the  domi¬ 
nant  frequencies  can  be  identified  with  inputs  that  are  not 
sufficiently  rich  to  excite  all  of  the  frequencies  (in  the  sense 
of  (GSl,LSl|).  Also,  a  few  interesting  issues  that  arise  due 
to  the  specific  nature  of  the  input  signal  can  be  discussed. 
The  frequency  of  the  input  sinewave  was  increased  gradu¬ 
ally  according  to  Figure  3  (this  graph  was  obtained  from 
a  simple  single  channel  lattice). 

All  four  sensors  see  these  frequencies  for  both  inputs. 


Figure  3:  Frequency  of  Input  Signal  as  Identified  by  Lattice 


Naturally,  the  noisy  sensor  (four)  gives  the  worst  estimates 
and,  in  one  exceptional  case,  can  not  identify  one  of  the  low 
modes.  In  all  of  our  studies  we  have  seen,  consistently,  that 
for  a  given  input  sequence  and  sensor,  high  order  filters  are 
needed  to  identify  the  frequency  while  good  prediction  can 
be  obtained  from  much  lower  order  lattices.  The  results 
below  are  all  obtained  from  a  lattice  filter  of  order  48  using 
the  sensor-two  measurements. 

Estimating  the  damping  ratios  is  another  topic  of  in¬ 
terest.  The  structure  is  widely  believed  to  have  non-linear, 
vibration-level  dependent  damping.  The  results  here  show 
that  although  a  linear  (and  proportional)  damping  model 
does  not  fit  the  data,  an  approximate  model  with  such 
properties  would  have  the  expected  damping  levels.  Also, 
it  was  observed  consistently  that  the  quare  wave  input 
(which  caused  substantial  vibration  in  the  grid)  results  in 
higher  damping  for  the  modes  of  the  structure. 

Three  sets  of  natural  frequencies  are  identified.  The 
first  set  is  a  group  of  three  low  frequencies  (including  the 
dominant  one  around  .56  Hz).  The  behavior  of  these  modes 
-the  third  one  in  particular-  is  quite  interesting.  The  sec¬ 
ond  and  third  modes  are  not  quite  identified  by  the  lattice 
at  the  beginning.  Once  these  modes  receive  substantial 
excitation  (i.e.,  when  the  input  frequency  is  close  to  the 
natural  frequencies),  the  lattice  results  in  good  estimates 
for  them.  This  can  be  seen  by  a  simple  comparison  of 
the  input  frequencies  (Figure  3)  and  the  natural  frequen¬ 
cies  (Figure  6).  Also  the  behavior  of  the  damping  ratio 
estimates  is  presented  for  one  of  the  modes  which  is  repre¬ 
sentative  of  the  rest. 

The  second  set  of  frequencies  are  those  which  the  higher 
modes  which  are  not  as  highly  excited  as  the  first  three 
frequencies  and  require  large  order  filters  (generally  speak¬ 
ing).  Figure  7  shows  the  first  three  of  these  modes. 

The  third  set  of  modes  are  a  group  of  high  frequency 
modes  which  have  extremely  low  damping  (typically  less 
than  one  percent  critical)  which  are  present  in  all  cases 
and  are  identified  very  quickly.  Further  investigation  of 
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still  data  (i.e.,  data  obtained  when  the  structure  is  not 
vibrating)  indicates  that  these  modes  are  from  the  sensors 
and  not  from  the  structure,  since  all  these  modes  were 
present  when  the  sensors  were  observing  the  still  grid.  The 
plots  in  Figure  8  show  the  behavior  of  three  of  these  modes. 

6.  Concluding  Remarks 

In  this  paper,  we  have  investigated  the  application  of 
lattice  filters  to  the  identification  of  flexible  structures. 
The  primary  advantage  of  lattices,  for  such  applications, 
is  the  flexibility  they  offer  with  respect  to  the  system  or¬ 
der.  By  being  order  recursive,  lattices  can  be  used  to  esti¬ 
mate  an  effective  order  for  the  system  and,  consequently, 
to  change  the  order  of  the  filter,  on-line,  without  any  loss 
of  data.  Experimental  data  from  the  NASA-Langley  grid 
was  used  to  show  some  of  the  results  obtained.  The  be¬ 
havior  of  the  prediction  for  different  order  filters  and  the 
corresponding  estimates  for  the  system  natural  frequen¬ 
cies  all  point  to  the  advantages  that  lattices  provide  for 
the  on-line  parameter  identifiction  of  flexible  structures. 

The  natural  frequency  estimates  shown  here  agree  with 
those  obtained  at  NASA-Langley,  see  [MWLN],  by  using 
extensive  off-line  techniques.  The  following  table  shows 
a  simple  comparisons  between  [MWLN]  and  the  results 
obtained  here;  i.e.,  [JGj. 


MWLN 

JG 

.3695 

**# 

.58 

.560 

1.046 

*•* 

1.420 

1.37 

2.128 

2.07 

The  lowest  mode,  at  .3695  Ha,  is  the  rigid  body  mode  due 
to  the  cables  connecting  the  grid  to  the  ceiling.  Since  this 
mode  was  not  excited  in  the  experiments  that  data  for  this 
paper  was  collected  from,  this  mode  is  not  identified  here. 
The  third  mode,  at  1.046  Hz,  has  nodes  at  the  sensor  loca¬ 
tions  that  were  used  here.  This  mode,  therefore,  could  not 
be  identified  with  the  sensors  used  for  the  data  presented 
here. 

The  accuracy  of  the  results  obtained  here  (less  that  5% 
variation  in  frequencies  shown  above,  for  example)  is  in¬ 
dicative  of  the  potential  new  applications  of  lattices  for 
on-line  identification  (and  adaptive  control  of)  of  flexi¬ 
ble  structures.  We  are  currently  investigating  using  white 
noise  as  the  input  (to  excite  all  modes)  and  tracking  of 
the  changes  in  the  parameters  when  the  mass  properties 
of  the  structure  are  altered  singnificantly.  Preliminary  re¬ 
sults  have  been  quite  encouraging.  Such  results  will  be 
presented  in  other  publications  in  the  near  future. 
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